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PREFACE 


This report consists of two parts. Part I is concerned with 
the static displacement and stress fields. Part II describes the 
dynamid . wheel-soil interaction. These studies were conducted 
during the period January 1, 1972 through October 31, 1972, under 
NASA Research Contract NAS8-25102 "Mathematical Characterization 
of Mechanical Behavior of Porous Frictional Granular Media," tech- 
nically monitored by Dr. N. C. Costes, The Geotechnical Laboratory 
of the Marshall Space Flight Center/ NASA, Huntsville, Alabama. 
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ABSTRACT 


A new definition of loading and unloading along the yield sur- 
face of Roscoe and Burland is introduced. This is achieved by noting 
that the strain-hardening parameter in the plastic potential function 
is deduced from the yield locus equation of Roscoe and Burland. The 
analytical results are compared with the experimental results for 
plate-bearing and cone-penetrometer problems and close agreements 
are demonstrated. 

The second part of the reports deals with the wheel-soil inter- 
action under dynamic loading. The rate-dependent plasticity or 
viscoelastoplastic behavior is considered. This is accomplished by the 
internal (hidden) variables associated with time-dependent viscous 
properties directly superimposed with inelastic behavior governed 
by the yield criteria of Roscoe and Burland. Effects of inertia and 
energy dissipation are properly accounted for. Exhaustive example 
problems are presented. 


vi . 
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PART I 


STATIC DEFORMATION AND STRESS FIELDS 


1-1. INTRODUCTION 


Recent achievements in the critical state soil mechanics 
advanced by Roscoe and others [1,2] have stimulated many other 
investigators searching for practical applications. Initial 
attempts have been made by Smith and Kay [3], Zienkiewicz [4], 

Chung and Lee [5], and Chung, Costes, and Lee [6] in the context 
of finite element techniques. The present study is an extension 
of [5,6] with some significant modifications in reference to 
interpretation of the yield criteria of Roscoe and Burland [1], 

In the previous works [5,6], the authors considered the 
strain- hardening parameter to be controlled by the constant yield 
stress, an independent material parameter, in addition to the basic 
material properties M,\and.K proposed by Roscoe and Burland [1]. 

However, in view of the fact that the equation of yield surface 
and subsequently the equation of yield locus as defined in [1] are 
based on the normality requirements of the plastic strain vector 

with strain-hardening phenomena ^incorporated in-the plastic poten — 

tial function, additional imposition of strain-hardening through a 
constant yield stress is unnecessary. Because the terms included 
in the plastic potential function [5,6] consists of deviatoric stress 
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invariant and the basic soil mechanics material properties (M,\, K ) 
associated with the mean pressure the later contributions in the 
plastic potential function must provide strain-hardening behavior in 
the sense of classical incremental theory of plasticity. This argu- 
ment leads to the standard manner of handling the plastic potential 
function in that the variation of the plastic potential function 
simply depends on the second deviatoric stress invariant and the 
strain-hardening parameter. If such variation is equal to zero we 
have a neutral loading, and the positive and negative values would 
indicate loading and unloading, respectively. The positive change of this 
potential function, therefore, shifts the yield locus in the devia- 
toric-mean stress space whose projection back to the void ratio - 
mean stress space lies entirely on the yield surface at all times. 

The constitutive relationships and the finite element equations 
are derived as demonstrated earlier [5,6], The plastic tangent stiff- 
ness matrix is updated for small increments of loading. The repeti- 
tive solution of the equilibrium equations continues until the total 
load is reached. Numerical examples for the plate-bearing and cone- 
penetrometer are presented to evaluate correctness of the procedure. 
Comparisons with test results indicate close agreements. 

1-2. YIELD CRITERIA AND PLASTIC STIFFNESS 

/ 

We record here the following basic assumptions of the criti- 
cal state theory: (1) the soil material is continuously distri- 

buted over its whole volume with its behavior described by a ma- 
croscopic model; (2) the mechanical behavior of cohesive and cohesion- 


less soil depends only on effective stresses independent of the presence or 
absence of pore pressures. The consequences of these assumptions lead to a 
complete description of soil behavior in a space of void ratio e, mean pres- 
sure p, and deviatoric pressure q. The deviatoric and volumetric strains 
corresponding to q and p along the yield locus are then related by means of 
the normality principle of plasticity theory as shovm in Figure 1. 

The mathematical model of pre- yield behavior may be based on the simple 
assumption of complete rigidity or elasticity, although some evidence exists 
of irrecoverable plastic shear distortion in this range [1]. For simplicity 
we may use the elasticity theory for the range of elastic wall (Figure 1) . 

To deal with irrecoverable volumetric and deviatoric strains and re- 
coverable volumetric strains we turn to the equation of yield locus, 

MS 


P_ 

Po 


= 0 


M a + "Hs 

where ”0 =q /p ; p 0 . is the mean pressure corresponding to q = 0; and M is 
the slope "0 at the critical state line, 

6 sin 9 


( 1 ) 


M = 


( 2 ) 


3-sin 9 

in which 9 is the angle of internal friction. 

The incremental plastic (irrecoverable) volumetric strain is 

C p ) 


dv 


( p ) 


de 


1 + e 

The overall void ratio change along the isotropic compression curve is 

d Po 


( 3 ) 


de = - X 


Po 


( 4 ) 


whereas the incremental recoverable void ratio is given by 
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de< r >=-H 

' Po- 

Here X and h are the compression index and swelling index, respectively. 
The incremental irrecoverable void ratio is then obtained from (4) and 
(5) as 

de< p) = - (X-H ) 

Po 


(5) 


( 6 ) 


At this point we introduce the equation of yield surface in the form 


[ 1 ] 


EL _ ( Mi-J \ 

p. # _ j 


C - *) 


(7) 


^M 3 + T) a- 

in which p e is the equivalent pressure corresponding to that void on the virgin 
isotropic consolidation line whose projection to the p - q space is P 0 • 
Therefore, setting p e = p 0 in (7) leads to 

a : (1 -h/x) 


Po 


X M 1 


( 8 ) 


Under triaxial compression, the second deviatoric stress invariant becomes 

J = 1/3 (a u -033 y = 1/3 q p (9) 

which gives 

q = 73T (10) 

Substituting (10) into (1) and rearranging yields 

3 J + p M s ( P - p 0 <) = 0 

or 3J - A s = 0 ' (11) 

where A 3 = •p'M s (!p - p 0 ) (12) 

It should be noted that (11) assumes the identical form as the plastic potential 
function F(J, A) in the sense of classical incremental theory of plasticity, 

F (J, A) = 3J - A 3 = 0 (13) 


$A 7^p 0 

The associated .flow rule for the incremental plastic volumetric strain 
dv * and the incremental plastic deviatoric strain tensor d:^,^ may be 


written. 

respectively. 




bA bp 

dX 

(17) 


,*(0 ftp bj 

dX 

(18) 

in which 

dX is the positive constant. Here dv l,} 

may also be expressed in 

an alternate form from (3) and (6), 




dv (P) 

1+8 

d:P 0 

(19) 

Equating 

(17) and (19) and using (16) give 
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dv 


( P ) _ ( X — K ) dp n 


(t? ' bF bA , \ 
_ -( x - h) IbJ dJ + bAbP dp J 


(Ue ). ..bib* (l + e) P hi hit hi hl^ 

U+e;p o bA b p U+e; P ° bA bp b A bp 0 . 


( 20 ) 


The incremental total plastic strain tensor is given by 

dV D ( r P> = di!;^ + 1/3 dvi n P) 6 on 


( 21 ) 


in which 6 mn i s the Kronecker delta. Using (17) through (20) in (21) yields 


where 


dV ( D P) = B mn R^ b .do 


jj * 


<*fB 


L bJ 1 fcFbA , \ M v 
( ba nn + 3 bA bp 6on j (X ' H) 


p (1+e) ^ 

P ° ; bA bp bA bp 0 


( 22 ) 


(23) 


also. 




hi 


hi. 


bJ dJ = 3 ^~ d « B " =^ d "° n 


dp = (2PM 5 - p 0 ’M ? ) dp 


bF bA 

^a^T = ‘ m p 


Substituting (25) into (23) gives 


Spn 3 ^ 


q n 


o n 


(24) 

(25a) 

(25b) 

(25c) 

(26) 


in which 


M 


= f- (2P - P 6 ) 


b = 3a(l+e)MVp 0 ^ /( X-h) 


rQ'P 


Similarly, R#pda p in (24) is given by 


(27) 

(28) 


where 


Su = 2 < 3 \\ - a sr -O 33 

Spp = 20pp -Oxx _p, 33 

S33 = 2a 33 -cjxx _ o 2 j> 

Sis = 60 j a , S S3 - 60 2 3 , S 3 I = 6o 31 

The incremental total strain tensor dV Bn is the sum of the incremental 

( p ) 
n ■*' 


elastic strain tensor dY and the incremental plastic strain' tensor dy n 
Therefore , 


dy^ <jv jy ( p ) 

u ' mu — uv tnn UT mn 


(30) 




The incremental total stress tensor dV is then given by 

CVp oPnn ...(e) /oi\ 

da a D ( e) dY^' (31) 

Q*P m n 


in which D( e ; is the standard elasticity matrix. Substituting (30) and 
(22) into (31) yields 

da"** * D ap '"(dY BB - B^Rjj da 1 i ) 

In view of (1.4b) and (24), and (32), we obtain 

K, . fn”” (dV„ - B n n do 1 ” 8 )] + |J ^ dP 0 . 0 


or 


from which 


Rn 


T.7 < dY -"- B - % 


Op 

do K ) 


i - R , 

J ctf 


do 


of 


(32) 


_ , OP 

R da = 

a P 

Substituting (33) into (32) gives 


where 


Rra P rBBn d Y 
1+ R rs B Bn D r “ 


orP 

f. CkP on 
D 

op 

an \ 

da = 

+ D 

dV n 


V <•) 

(p 

> J ° 

O' P on 

k r 

O' p n 

D 

Bki R 

ii D 1Jon 

IP) 

1 + B 

re ®tu 

q r#tu 


(33) 

(34) 

(35) 
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which is identical to the form obtained by the authors earlier [5,6]. 
Now, the yield criterion equation (14) is written as 

dF = R - M% dp 0 (36) 

where dp 0 can be determined from (8), 


dp 0 = g dp + h S t J d aiJ 


(37) 


in which 




(-H/x) (38a) 


h *.«• i 1 ' f)( 1 + fp^) 


(38b) 


Substituting these in (36) yields 


dF 


-[ s . 


Q + a 6 _ - M 
p aft 


(5 « V ♦ 


( 39 ) 


which is then used for determining the status of loading, neutral 
loading, and unloading as defined in (15 a, b, c). 
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1-3. APPLICATIONS 

1.3,1 Plate Bearing 

Based on the definition of yielding given by (15) the finite 
element computer program was written to solve boundary value prob- 
lems. The program listing and data input format are given in 
Appendix 1. and Appendix 2 , respectively. 

Figure 2 shows the geometry of a plate bearing problem. The 
load-displacement curves for center of plate are shown in Figure 3 
comparing the experimental results of Namiq [8], It should be noted 
that the plane strain conditions of Namiq's experiments with a square 
box are approximated here in the analysis by an equivalent axl- 
symmetric cylindrical box. The material constants given by Namiq 
are angle of internal friction = 35°, initial void ratio ® = 0.875, 
initial density V = 0.0147 N/cc. Other constants needed in this 
analysis are listed in Figure 3. It is seen that the load-displacement 
curve for the compression index \ = 0.05 follows very closely the ex- 
perimental results whereas \ = 0.13 gives slightly larger displace- 
ments. It is interesting to note that from the void ratio-pressure 
curves given by Namiq the compression index can be estimated indeed 
to be approximately 0.05. Here the swelling index h = 0.003 is used 
for both cases. For elastic behavior the soil modulus E g = 10 N/cm^ 
and Poisson's ratio v 8 = 


.45 are used. 
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Axis of Symmetry 



Figure 2: Plate Bearing Geometry 

» ■ , 
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Deformed shapes of the finite elements adjacent to the bearing 
plate are shown in Figure 4 for the loading increments at F = 2.5 N/cm® 
and F = 5 N/cm s . These results correspond to \ * 0.05 which gives 
the same displacement at the center of plate'as Namiq. Unfortunately, how- 
ever, no further comparison can be made as Namiq does not show such deformed 
shapes in his experimental results. 

1.3.2 Cone-Penetrometer 

The geometry for a cone-penetrometer problem is shown in Figure 
5. Experiments for the cone -penetrometer were undertaken and the test 
set-up' is shown in Figure 6. Both smooth and rough aluminum cones 
were used and loaded through the lunar soil simulants under the strain- 
controlled loading devices. These measurements are plotted in Figure 
7 and compared with analytical results. The axisymmetric interface 
elements developed by Chung and Lee [5] are used to model contact 
areas between the cone and soil. Because of the lunar soil simulants 
being extremely soft compared with the metal cone the shear modulus 
and rotational modulus for the interface elements were set equal to 
zero. Experimentally determined material constants for the lunar 
soil simulants used in the tests are also given in Figure 7. The same 
material constants were used in the analysis with the exceptions of 
soil modulus E 0 = 10 N/cm® and Poisson's ration y,. = .45. The anal- 
ytical solution gives results somewhere between the rough and smooth 
cones. 

The deformed geometry of soil is shown in Figure 8. For exces- 
sive alterations of finite elements in shape it would appear that 


5N/cm a 
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renumbering of nodes is necessary to update the stiffness matrix 
based on new geometry. It is believed that such treatment would 
improve the solution considerably. 

1-4. CONCLUSIONS 

A new definition of loading and unloading along the yield sur- 
face of Roscoe and Burland is introduced. This is done by noting 
the strain-hardening parameter in the plastic potential function. 

With the differential of the plastic potential function with respect 
to the second deviatoric stress invariant and the strain-hardening 
paramenter being positive or negative the manner of loading and un- 
loading is clearly determined. This is an improvement from the pre- 
vious definition of yielding through a constant yield stress. 

The forms of plastic stiffness matrix and the finite element 
equations, however, are unchanged. Applications of the present anal- 
ytical formulation to a number of boundary value problems are pre- 
sented. The analytical results for the plate bearing and cone- 
penetrometer problems indicate good agreements with the experimental 
results . 

Our ultimate goal is to characterize the material parameters of 
the lunar soil. Such a task depends on correct constitutive relation- 
ships and a computational scheme which provides the results of load- 
deformation. With this facility available exhaustive computer runs 
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for various combinations of material constants are to be compared 
with the data brought back from the lunar exploration. To this end 
the present study has provided the basic analytical tool to prepare 
for such an undertaking. 
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COMPUTER PROGRAM LISTING 
(Static Analyal® - A*lsy*trlc) 



liELIfluL UtCK t »SYS'MJP 
ELI 005-11/24-14:39 

000001 UOO l»ELT » SIH NASA + TPFS. . MA IN . . , 1326561330 1 0 

0OUUO2 UOU L OOOOOIOO 

UUOU03 UOU C 00000200 

0UUUU4 UUU C THE FINITE ELEMENT ANAL1SYS OF AXISYMMETRIC SOIL MEDIUM 00000300 

0U0005 UOO C BY A SOIL PLASTICITY THEORY 00000400 

UUU0U6 UUU C 00000500 

UUU007 UOU 00000600 

000008. UOU C 00000700 

000009 ,000 PARAMETER NODS=30u»NELS=260.NF=20000*MAX=600 00000800 

UOUOXO UOU C 00000900 

uuuoil UOO COMMON /bLKO/ T I TlE ( 20 ) « INODE , NELEM > NAPC , NBC » NINCR , NC YCL tEPSLON 00001000 

UUUU12 UOU COMMON /6LKJ/ W ( 6 ) . H ( 6 ) . AR ( 4 ) . BR ( 4 > . CR ( 4 ) , A? ( 4 ) , BZ ( 4 ) . CZ ( U ) • 00001100 

0UUU13 UOU * BN(4) »CN<4> »DN<4 J . TYPEA (4« 4 )". TYPEB (4 »4 > rTYPEC <4»4> *TYPEE(4»4 ) , 00001200 

00U014 UOU * TYPEF(4»4) *TYPEtiI4»4) fAO»bO»CO.RT.RB»RA»RCf IC»JC,KC,LC.NEL 00001300 

UUUOlb UOU COMMON /BLK2/ ID ( NODS . 2 > . I JKL < NELS • 4 ) ,DE1 , DE2 00001400 

UUU016 UOU COMMON /tjLK3/ XK ( Np ) . APF ( MAX ) t IMAX » IHP . IhB I » LT t LAST 00001500 

U0U017 000 COMMON /BLK4/ STRS(NELS#4) » DM AT (NELS » 4 » 4 ) .DELMM). (MAX) »POP 00001600 

OOUOIO UOU COMMON /bLK5/ LIE ( NELS # 4 > 4 ) *SlG8A(N£LS) » OS 1 66 A ( NELS ) >DELL>YSTRSf 00001700 

0UU019 UOU 1 FInC.FN.ULOAD»FEl*PMAX»L)LMAX 00001800 

00U020 UOO COMMON /bLK6/ SIGK # SIGZ * S IGT * TAUZR *0 (4 » 4 ) rSTIFF(6»8) *KK(NELS»8) , 00001900 

O0U021 UOU 1 RINUOS) iZLNODSl . iOTUIS(MAX) 00002000 

00U022 OOU COMMON /dLK 7/ LiSTKS < NELS > u X. ARM < NELS» 4 > » AZM < NELS > 4 > » RTT < NELS > » 00002100 

000023 000 1 AOJINELS) 00002200 

00U024 UOO COMMON/ BLKO/ 1 NCR »PUEPTH(nELS) » VOIDI » AL AMDA ( DEPTH * PP ! 00002300 

0U0025 000 COMMON /bLK9/ P I . SMALLK . CK * BETA , PO > NFREE » NELST » ICASE , NR IGn > 00002400 

00U026 000 COMMON /bLKlU/ FRCM20*8. 8) ,IR(20*8»8) »XXL(20) »DZZ(20) »DRR(20) » ’ 00002500 

UUUU27 OOU 1 POR ( 20 ) » OELTS • EC * XNLC t DELD I ! 00'002600 

00UU28 UOU • COMMON /bLKll/ VOID ( NELS ) .OGAM ( NELS . 4 ) 00002700 

000029 UOU COMMON /BLM2/ SIUMX (NELS ) #OEP (NELS ) .EP ( NELS ) ► DEQST2 (NELS) »ES 00002800 

00UU30 UOU COMMON /GMTRY/ RO(NODS) »ZO<NODS> 00002900 

00U031 UOU C 00003000 

00U032 UOO C 00003100 

UU0U33 UOO NT APE = 2 00003200 

0UU034 UOU C 00003300 

0UU035 UOU C • 00003400 

O0U036 OOU C 00003500 

000037 UOO CALL SETUP 00003600 

00003b UOU C 00003700 

000039 000 L .INITIALIZES NECESSARY CONSTANTS FOR INTEGRATION SCHEME. 00003800 

U0U040 UOO C 00003900 

000041 UOU C 00004000 

000042 UOO CALL INPUT T NT APE ) 00004100 

000043 000 ISHEAR =1 00004200- 

000044 000 . C 00004300 

000045 UOO C 00004400 

000046 000 FN = FEL 00004500 

000047 OOU C 00004600 

000048 UOO — 00004700 

000049 000 C 00004800 

000050 UOU C START MAIN ITERATION- LOOP. 00004900 

000051 UOO C 00005000 

000052 UOU C . — — 00005100 

0U0053 UOU DO 990 NI = 1. NInCR 00005200 

000054 .000 C 00005300 

UUUU55 UOO FN = FN + FINC 00005400 


UUUUS6 

000 


iter = o 



00005500 

U0OU57 

OOU 

c 




00005600 

00UU56 

000 


IF(Nl.EQ.l) GO TO 950 



00005700 

UUUUS9 

000 


IShEAR = Ni 



00005800 

UUUU60 

OOo 


CALL Z£R0(XK.NF»1) 



00005900 

UUU061 

000 


REtvIND NT APE 



00006000 

U00U62 

000 

c 




00006100 

0UUU63 

000 


00 940 NlL = 1 t NELEM 



00006200 

UUU064 

000 


IL = I JKL (NEL » 1 ) 



00006300 

UUUU65 

000 


Jt = 1 JKL ( NEL » 2 ) 



00006400 

0UUU66 

000 


KL = 1 JKL ( NEL » 3 ) 



00006500 

U0UU67 

000 


Lo = 1 JKL (MEL * 4 ) 



00006600 

U 0 U 0 6 8 

000 


NF 1 = NEL - NRiGU 



00006700 

UUUU69 

000 


IF(NEL.LE.NRIGU) GO TO 938 



00006800 

UOU07U 

000 


SIGZ = SrRS(NEL.l) + DSTRS<NEL»1) / 2. 



00006900 

0UU071 

000 


SIGN = S 1 Hb ( nEl » 2 ) + 0bTKS(NEL»2> / 2. 



00007000 

(JUU072 

000 


SIGl = STRS(NEL.j) + ObTRS ( NEL * 3 ) / 2. 



00007100 

UO0U73 

000 


TAUZR= S 1 Rb ( NEL ? 4 ( + ObTHS ( NEL * 4 ) / 2. 



00007200 

0UU074 

000 


IF(NbC.NE.O.ANU.Nl-T.LE.NtjC) GO TO 699 



O0OQ730O 
on n n74 n n 

00UO76 

000 

c 




00007500 

UUU077 

000 


CALL DMATRX(NI) 



00007600 

0UU076 

uoo 

c 




00007700 

U0U079 

000 

c 

CALCULATE STRESS DEPENDENT MATERIAL 

PROPERTY 

matrix (D>. 

00007800 








000081 

000 

c 




00008000 

0OU082 

000 


GO TO 837 



00008100 

000083 

ooo 

c 




00008200 

000084 

000 


938 CO 400 1 = 1»4 



00008300 

00008b 

000 


DO 4U0 J = 1 »4 



00008400 

00008b 

ooo 


OMAT ( NEL » I r J ) = UE(MELrlfJ) 



00008500 

000087 

ooo 


400 CiliJ) = DMA! (NEL# if J) 



00008600 

000088 

ooo 


837 CALL ST IF F2 ( Nl »NT«PE ) 



00008700 

00U089 

ooo 


,GC TO 838 



00008800 

000090 

ooo 

c 




00008900 

000091 

ooo 


e99 CALL FRICTN(lC#KCfNEL# NFT # 1 SHEAR > VO ID I ) 



00009000 

000092 

ooo 

c 




00009100 

000093 

ooo 


838 CALL ASSEMb(NEL.NFT) 



00009200 

000094 

ooo 


940 CONTINUE 



00009300 

00009S 

ooo 


950 CONTINUE 



00009400 

000098 

ooo 


ITER = ITER + 1 



00009500 

00O097 

ooo 


I NCR ITER 



00009600 

000098 

ooo 

c 




000097,00 

000099 

ooo 


CALL DISRL ( NFREE » Ni , 1N00E dNCR ) 



00009800 

000100 

ooo 

c 




00009900 

000101 

ooo 


00 340 1 = l.INOUt 



00010000 

000102 

ooo 


JJ = I * 2 



00010100 

000103 

ooo 


II = JJ - 1 



00010200 

000104 

ooo 


Z(l) = ZOd) + TOlUlS(lI) + XK < L AS 1 + 11) 



00010300 

00010b 

ooo 

340 R(l> = Rud> + TO t DIS ( JJ) + XMLAST+JJ) 



00010400 

000108 

oou 

c 




00010500 

000107 

OOU 


CALL STRAIN ( NI * ITcR ) 



00010600 

00O108 

oou 

c 




00010700 

000109 

ooo 

c 

Summing OF STRLbSES and displacements 

FOP EACH 

INCREMENTAL 

STEP00010800 

00O110 

000 


DO 310 ITT = l.NFhEE 



00010900 

000111 

ooo 


310 TOTOISdT T ) = TOluISdTT) + XK ( I TT+LAST ) 



00011000 

000112 

uoo 


DO 329 I = 1 f NELEm 



00011100 


UUU113 
oouiit) 
OUUUb 
0U0116 
U0U117 
OOU118 
U0U119 
0 OU 120 
00U121 
000122 
000123 
000124 
000125 
000126 
000127 
000128 
000129 
000130 
000131 
000132 
000133 
000134 
000135 
000136 
000137 
000138 
000139 
000140. 
000141 
000142 
000143 
000144 
000145 
000146 
000147 
000148 
000149 
000150 
000151 
000152 
000153 

000154 

000155 

000156 

000157 

000158 

000159 

000160 

000161 

000162 

000163 

000164 

000165 

000166 

000167 

000168 

000169 


000 


DO 329 J = 1 #4 

00011200 

oou 

329 

STKS(I#J) = STRS ( 1 » J ) + DSTRS(I.J) 

00011300 

000 

c 


00011400 

000 


WR 1 TE ( 6 . 689 ) NlNCl<»Nl»FN 

00011500 

000 


CO 320 I = l.INOOe. 

00011600 

oou 


11 = (1-1) *2+1 

00011700 

000 


Jj = 11 + 1 

00011800 

oou 

320 

WR I TE ( 6 # 690 ) I # ( TOTO IS (LX ) » LX = II#JJ) 

00011900 

000 


WRlTt(6» 691 ) 

00012000 

OOU 


DO 330 I = 1#NELEM 

00012100 

uoo 

330 

WRITE (6,692) I# ( S IKS ( I » U) #U=1#4) 

00012200 


000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 


990 


500 

600 

689 

690 

691 

692 

693 

694 

695 


UELT. 

C 

C 


IFIFIm.GT.PMAX) STOP LOADMX 00012300 

CONTINUE 00012400 

00012500 

00012600 

FORMAT (1015) 00012700 

FORMAT ( // * ITER ' # 15 » ' OMAX ' # E 12 . 5 » ' 0E2'#E12.5> 00012800 

FORMAT (1H1#10X«' TOTAL DISPLACEMENT NO. OF INCREMENTAL STEPS 00012900 

*=' #2I5//5X# 'NODE* #5X» *2 - OISPL • » 20X# 'R - DISPL* #5X» *FN = • #E12. 5/ ) 00013000 
FORMAT ( 4X »I5»E15.7#10X#E15.7) 00013100 

FORMAT (///#10X» * TOTAL STRESSES • //5X #* ELEM •# 5X #• SIGMA - Z*,T28# 00013200 

l'SIGMA - R' ,T42. • TANGENTIAL* »T58# 'TAU - ZR') 00013300 

FORMAT(Id» 4F14 .6 ) 00013400 

FORMAT ( // ' TAU# SIG. RAT# DELTS# ISHEAR# • .4E12.5# 15//) 00013500 

FORMAT < //2UX # 'NEw GEOMETRY AT THE END OF LOAD INCR.'#I5//) 00013600 

F OHM AT(I10»2P 15.6) 00013700 

STOP 00013800 

END 00013900 

SIH NASA*TPFS. INPUT# . #132661133010 
SUBROUTINE INPUT (NTAPE) 


PARAMETER N00S=3UU # NELS=260 . NF=20000 # M AX=600 


COMMON /BLKO/ T I TlE ( 20 ) # INODE . NELEM »N APC , NBC »NINCR # NCYCL #EP5L0N 
COMMON /BLK1/ W ( 6 J » H ( 6 ) # AR < 4 ) #BR ( 4 ) » CR ( 4 ) , AZ ( 4 ) , BZ ( 4 ) » CZ ( 4 ) » 

* BN (4) #CN(4) »0N(4) #TYPtAt4#4> »TYPEB(4#4) #TYPEC(4#4) »TYPEE(4.4> > 

* TYPEF 1 4 » 4 ) . T YPEG 14 » 4 ) # AO # bO # CO ,RT # RB # RA » RC # IC # JC #KC ,LC # NEL 
COMMON /BLK2/ 10 ( NODS # 2 ) # IUKL (NELS # 4 ) #DE1 # DE2 

COMMON /8LK3/ XK (Np ) » APF l MAX ) # I MAX » I HR # IHBI # LT # LAST 
COMMON /fc>LK4/ STRS (NELS # 4 ) • DMAT (NELS # 4 # 4 ) #DELNM1 (MAX ) # POP 
COMMON /BLK5/ UE (NELS#4.4) » SlGB A (NELS ) #DSIGBA ( NELS ) #DELL»YSTRS» 

1 F1nC#FN.UL0AD#FE>_»PMAX.ULMAX 

COMMON 7bLK6/ SIGN #SIGZ #SIGT » TAUZR #D (4#4) »STIFF (8.8) #KK(NELS#8) # 
1 R ( NODS ) » Z ( NODS ) # I OTD IS ( M AX ) 

C0MM0N/BLK8/ INCR . POEPTH ( NELS ) > VOIDI # ALAMDA # DEPTH # XMS 
COMMON /BLK9/ PI # SMALLK » CK »BET A #P0 #NFREE »NELST # IC ASE #NRI6D 
COMMON /BLK10/ FROM 20 # 8 # 8 ) . TR ( 20 # 8 # 8 ) » XXL ( 20 > #DZZ (20 > # DRR ( 20 ) » 
l POR (20 ) » DELTS# EC » XNUC #DELD 
COMMON /BLK11/ VOID ( NELS) »OGAM (NELS #4 ) 

COMMON /bLK12/ SIGMX (NELS )# DEP (NELS > #EP (NELS) # DEGST2 (NELS ) ,ES 
COMMON /GMTRY/ RO (NODS )# ZO (NODS ) 

REWIND NTAPE . 

READ ( 5 > 510 ) (TITLE(I)»I=1.20) 

READ (5# 500) INODE »NELEM »NAPC * NBC #NINCR# NCYCL » ICASE»NRIGD»NULOAD 
REAO(5#530) YSTRS#DELL# ZETA#RMAX #DI MAX 
READ ( 5 • 511 ) DZI.ORl 


00000100 

■00000200 

00000300 

■00000400 

00000500 

00000600 

00000700 

00000800 

00000900 

00001000 

00001100 

00001200 

00001300 

00001400 

00001500 

00001600 

00001700 

00001800 

00001900 

00002000 

00002100 

00002200 

■00002300 

00002400 

00002500 

00002600 

00002700 

00002800 


000170 

uuu 

L 


021 = SHEAR MOD. FOR INTERFACE ELEMENTS. 

00002900 

000171 

0(JG 

c 


OKI = NORMAL FOR INTERFACE ELEMENTS. - 

00003000 

00017? 

000 



REAO (5*511) EC.XNUC.ES.XNUS 

00003100 

000173 

000 



READ (5. 530) PI.SMalLK.XI.VOIOI.PO.DEPTH. ALAMDA.EPSLON 

00003200 

000170 

000 



WR1TE(6*S3U)PI* SMALLK *Xl*VOIDI*PO* DEPTH * ALAMO A * EPSLON 

00003300 

00017b 

uoo 



WRITE (6*630 ) YSTRbrDELL 

00003400 

000176 

000 



WRITE (6*631 ) ZET A * 0? I > DR1 

00003500 

000177 

ooo 


631 

FORMATt//* ZE TA* MOD. FOR INTERFACE ELEM.* SHEAR * NORMAL* 

00003600 

000176 

000 


] 

L 3F15.6) 

00003700 

000179 

000 



POP = EXP 1 1 . -SMALi-K/ ALAMOA ) 

00003800 

oooleo 

Ooo 



II = 0 

00003900 

oooiei 

ooo 



00 101 I = 1 * INOOE 

00004000 

000182 

ooo 



00 101 J = 1*2 

00004100 

000183 

ooo 



II = II + 1 

00004200 

000184 

ooo 


101 

ID ( I * J) = II 

00004300 

000186 

uoo 



BETA=ALAMOA-SMALLK 

00004400 

000186 

ooo 



PI=PI*3. 14159/180. 

00004500 

00U187 

ooo 



EPSLON = EPSLON * 3.14159 / 180. 

00004600 

000188 

ooo 



CELTS = TAMEPSLON) 

0000,4700 

000189 

ooo 



PI=SIN(PI> 

00004800 

000190 

ooo 

c 


PI IS SINE (PI) 

00004900 

000191 

ooo 



XM = 6 . *PI/( 3 .“Pi ) 

00005000 

000192 

ooo 



XMS = XM * XM 

00005100 

000193 

Ooo 



WR 1TE ( 6 * 60u ) (1 ITLt(I) .1 = 1*20) 

00005200 

000194 

ooo 



00 lUO I = l.lNOOt 

00005300 

000195 

uoo 



RE AO ( 5 * 52 0 ) Z ( 1 > * K ( I) . IZ . IR 

00005400 

000196 

ooo 



ZO(I) = ZU> 

00005500 

000197 

ooo 



RO ( I ) = HU) 

00005600 

000198 

ooo 



IFU2.NE.0) 10(1*1) = 0 

00005700 

000199 

uoo 



IF(Ik.ME.O) IDU.2) = 0 

00005800 

000200 

uoo 


100 

WR ITE ( 6 * 620 ) IiZUliHIl) rlZ.IR 

00005900 

000201 

000 



WR1TE<6«50I) 1N0UE * NELEM * NAPC * N1N0R * MC YCL 

00006000 

000202 

ooo 



NFREE = INODE * 2 

00006100 

000203 

uoo 

L 



00006200 

000204 

ooo 



WR ITfc. ( 6 * 651 ) 

00006300 

000205 

ooo 



REAO (.6*640) ( ( IJKtUNEL* J) *0=1*4 ) *NfcL=l .NELEM) 

00006400 

000206 

ooo 



WRITE (6. 650) (NEL* UUKL<NEL*J> *0=1*4) * NEL=1 * NELEM ) 

00006500 

000207 

ooo 

c 



00006600 

000208 

Ooo 

c 


FIND HALF BAND WiUlH AND ACTUAL SIZE OF MATRIX (XK) 

00006700 

000209 

ooo 

c 



00006800 

000210 

ooo 



IMAX. = 0 

00006900 

000211 

ooo 



00 81)0 NEL = 1 * NEi_EM 

00007000 

000212 

Ooo 



00 700 1 = 1*4 

00007100 

000213 

ooo 



IN = I JKL (NEL* l ) 

00007200 

000214 

ooo 



KK(NtL.I) = ID ( IN * 1 ) 

00007300 

000215 . 

ooo 


700 

KK(NEL.I+4) = 1DUN.2) 

00007400 

000216 

ooo 



00 7999 I = 1*2 

00007500 

000217 

oou 



11 = J + 1 

00007600 

000216 

ooo 



00 7999 J = II. 4 

00007700 

000219 

ooo 



IDIF = IOKL(NEL.I) - IOKL(NEL.J) 

00007800 

000220 

ooo 



IFUDIF.LT. 0) iniF = -IDIF . 

00007900 

000221 

ooo 

7999 

IFIIUJF.ST.IMAX) IMAX = iOIF 

00008000 

000222 

POO 


800 

CONTINUE 

00008100 

000223 

ooo 



IF(NBC.NE.U) RE AD ( 5 * SOU ) ( UO ( I * 0 ) * J=1 , 2 ) . 1=1 * NRC ) 

00008200 

000224 

ooo 

c 


IMAX = MAX DIFFERENCE IN AOJACENT NODE NO. 

00008300 

000225 

ooo 



IHb = UMAX + 1) * 2 

00008400 

000226 

ooo 



IHDI = IHB - 1 

00008500 



IJVU227 

000 


LT = IHS * IHB1 / 2 


00008600 

000226 

000 


LAST = LT + (NFREt - IHBI) * IHB 


00008700 

00U229 

000 


WR I T£ ( 6 * 640 ) IMAX * 1H6 * LT * L AST 


00008800 

UUU231 

000 


IF(NRIGD.NE.O) CAi-L ELASTC < D* NR IGD • EC . XNUC > 


00009000 







000233 

000 


NLST = NELEiy - NR1G0 - NbC 


00009200 

000234 

oou 


CO 900 NEL = 1 * NEt£M 


00009300 

00U235 

000 


1C = IUKL(NEL.l) 


00009400 

0OU236 

000 


JC = I JKL (NEL#2 ) 


00009500 

000237 

000 


KC = I JKl ( NEL # 3 ) 


00009600 

000238 

000 


LC = I JKL f NEL * 4 ) 


00009700 

0UO23S 

000 


IF ( NEL • Lt . NR IGu ) GO TO 696 


00009800 

000240 

000 


NFl = NEL - NRIGU 


00009900 

000241 

000 


IF(NKT.GT.NPC) GO 10 896 


00010000 

000242 

000 


D2/<nFT> = D2I 


00010100 

000243 

000 


DRh(NFT) = OKI 


00010200 

000244 

000 


CALL FRICTN( IC * KC *NEL * NFT . 0 » VOIDI ) 


00010300 

000245 

000 


GO TO 899 


00010400 

00024b 

000 


896 IF ( ILASE • NE. U) GO TO 897 


00010500 

000247 

000 

c 



00010600 

000248 

000 


CALL AREAA(IC.JC,nC,LC*AREA> 


00010700 

000249 

ooo 


OOIOR = VOlOl 


00010800 

000250 

000 


VOIC(MEL) = V010R 


00010900 

000251 

000 


DELNR1 ( NEL ) = AREA 


00011000 

000252 

000 

c 



00011100 

000253 

000 


897 NElSI = NR1G0 + NbC + 1 


00011200 

000254 

oou 


IFlNEL.EO.NELSl )CALL ELASTC ID rNLSI >ES » XNUS) 


00011300 

000255 

000 

c 



00011400 

000256 

ooo 


896 DO 111 I = 1*4 


00011500 

000257 

000 


DO 111 J = 1*4 


00011600 

000256 

ooo 


DE ( NEl. » I * J ) = U(I*J> 


00011700 

000259 

ooo 


111 DMAT (NEL * I *U) = Oll.J) 


00011800 

000260 

ooo 

c 



00011900 

000261 

ooo 


CALL STIEFKnTAPE) 


00012000 

000262 

ooo 


CALL STIFF2 ( 0 * NT AFt ) 


00012100 

000263 

ooo 


899 CALL ASSEMd ( NEL * NET ) 


00012200 

0GU264 

ooo 

c 



00012300 

000266 

ooo 


900 continue 


00012500 

000267 

ooo 


IF(NAPC.nE.O) call PTLOAQ(NAPC*ULOAD) 


00012600 

0U0266 

OOO 


901 IF(NULOAO.NE.O) Call EOLOAD(ULOAD) 


00012700 

000269 

uoo 

c 

« • 


00012800 

000270 

oou 


SCAL = N1NCR 


00012900 

000271 

ooo 


DO 200 I = 1. NFREt . 


00013000 

000272 

ooo 


200 APE ( 1 ) = APF 11) / SCAL 


00013100 

000273 

uoo 


FEL = 0. 


00013800 

000274 

ooo 


FINC = ULOAD / SCAL 


00013300 

000275 

ooo 

c 



00013400 

000276 

ooo 


500 F ORMAT ( 10 1 S ) 


00013500 

000277 

oou 


501 FORMAT*//' NUMBER OF NODES =».I5/» NUMBER OF 

ELEMENTS = •* 

100013600 

000276 

ooo 


15/ • NUMBER OF APLIED CONCENTRATED LOADS =»*I5/» 

number of 

0001^700 

000279 

ooo 


2INCREMENTAL LOAD STEPS ='.I5/* NUMBER OF ITERATIONS PER EACH 

100013800 

000280 

ooo 


3NCREKENTAL LOADING =»*I5/) 


00013900 

000261 

uoo 


510 FORMAT (2UA4) 


00014000 

000282 

ooo 


511 FORMAT (4F20.5) 


00014100 

000263 

ooo 


520 FORMAT (2F 10, 4*215) 


00014200 


00U284 

0UU28S 

OUU286 

000287 

OU0286 

OUU289 

00U29R 

000291 

0UU292 

000293 

00029*4 

000295 

000296 

000297 

000298 

000299 

000300 

000301 

000302 

000303 

000304 

000305 

00O306 

000307 

000308 

000309 

000310 

000311 

000312 

000313 

000314 

000315 

000316 

000317 

000318 

000319 

000320 

000321 

000322 

000323 

000324 

000325 

000326 

000327 

000328 

000329 

000330 

000331 

000332 

000333 

000334 

000335 

000336 

000337 

000338 

00033-9 

000340 


000 

000 

000 

ono 

ouo 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

uoo 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

oou 

000 

000 

000 

000 

000 

000 

000 

000 

000 


530 FOKMAT (8810.4) 00014300 

531 FORMAT!/' SINE. PHI ='*F1U.4*' SMALLK ='.F10.4*» KI ='*00014400 

*F10. 4* ' POROuSITY =' ,F10.4/» BET A = • * FI 0 .4 * • OVERBURD00014500 

♦EM PRESSURE ='*F1U.4*» DESIRED ACCURACY ='*F10.4*» PERCENT. '/00014600 

* * DtPTH OF SoIL MEDIA ='.F10.4/) 00014700 

540 FORMAT (415) 00014800 

600 FORMAT (1h 1*20X*20A4///30X* 'COORDINATE VALUES’ //Til * ’ NODE' *T30. '2-C00014900 

♦OORO* ,T5u> 'R-COCR^' ,T65. 'O.IF FREE TO Z'»3X*'0*IF FREE TO R«//> 00015000 

620 FORMAT(T10*Ib*I25>F10.4*T45.F10.4*T65*2(l8*5X)) 00015100 

630 FORMAT ( // • YltLD STRESS ='*E12.5*» DELL ='*El2.5//> 00015200 

640 FORMAT!//' I MAX * 1HB . lT * LAST » *4 1 10// ) 00015300 

650 FORMAT (517) 00015400 

651 F0RMATUH1* 10X * ' CnNECTI VITY • / ) 00015500 

RETURN 00015600 

END 00015700 

WELT *SIH NASA^IPFS.DMATkX* ** 132671133010 

SUBROUTINE DMATPX(Nl) OOOOOIOO 

t_ — — — 00000200 

PARAMETER NOUS=3OU*N£LS=26O.NF=2OOOO*MAX=6O0 00000300 

■ 00000400 

COMMON /dLKO/ T ITLE ( 20 ) • INODE * NELEM . NAPC > NBC * NInCR * IPRINT * EPSLON 00000500 
common /oLKi/ w(6i* h(6) *ar( 4) .fiRU) » cr(4) ,az(4) *bz( 4) *cz<4) * 00000600 

* Bn(4) *CN(4) * DN ( 4 ) *TYPEA(4*4) * T YPEB (4*4) * TYPEC ( 4*4) *TYPEE(4*4) * 00000700 

* TYPEF(4*4) *TYPEG (4 *4) > AO»bO*CO*RT *RB*RA*RC* IC * JC*KC,LC *NEL 00000800 

common /blk4/ sths(nels*4) *dmat<nels*4*4) *delnmi (max) .pop 00000900 

COMMON /8LK5/ BE ( NELS . 4 * 4 ) * SIGB A (NELS ) * DSIGBA ( NELS ) * DELL ► YSTRS * l 00001000 

1 FINC.FN.ULOAD.FEL.PMAX.DLMAX ! 00001100 

COMMON /dLK.6/ SI6R * SIGZ * S I6T * TAUZR *D ( 4 * 4 ) , STIFF ( 8 * 8 > * KK (NELS*8 ) * 100001200 

1 R(MODS) .Z(NOOS) * lOTOIS(MAX) 100001300 

COMMON /6LK7/ USTRS (NELS • 4 ) .BUM (NELS* 10 ) I 00001400 

COMM0M/BLK8/INCR.HuEPTH(NELS) *V0IDI.ALAMDA,DEPTH*XMS • 00001500 

COMMON /BLK9/ PI* SMALLK *CK * RET A *P0 *NFREE * NELST * IC ASE *NRIGo 00001600 

COMMON /BLK1 1/ VOlO(NELS) *DGAM (NELS * 4 ) 00001700 

COMMON/BLK 12/ SIGmX (NELS ) * DEP( NELS) *EP ( NELS ) * DE0ST2 (NELS) *ES 00001800 

DIMENSION RLR(4) *UB(4) *R0(4) 00001900 

00002000 

V01DR = VOID (NEL ) 00002100 

P = (SIG2 + SIGT + SIGR) / 3. 00002200 

TJ = ( (SlGZ-SIGR)^42+(SIGR-SlGT>^2+(SIGT-SlGZ)*^2)/6.+TAlJZR**2 00002300 
PSO =. P^P 00002400 

ETS =’ 3.*TJ/PSU 00002500 

POW = l.-SMALLK/ALAMOA 00002600 

SIGMX(NEL) = P» ( ( XMS+ETS ) /XMS) ** POW 00002700 

POP = SIGMX (NEL) 00002800 

AA = XMS ♦ ( 2 . ♦P-POP ) /.. 3. 00002900 

B6 = 3. ♦ AA ♦ XMS ♦ P ♦ POP * (l.+VOIDR) / BETA 00003000 

CO 100 I = 1*4 00003100 

DO 100 U = 1*4 00003200 

100 0(1*0) = 0E(NEL*I*U) 00003300 

SZ2 = 2.^SIGZ-SIGR-SIGT 00003400 

SRR = 2.^SI0R-SIGZ-SIGT 0001)3500 

STT = 2.^SIGT-SIG2-SIGK 00003600 

SZR = 6. ♦TAUZR 00003700 

RLB(l) = SZZ + AA 00003800 

RLb (2 ) = SRR + AA 00003900 

RLB(3) = STT + AA- 00004000 

RLb (4 ) = SZR 00004100 


000341 

000 


CELP = (uSTRS(NEl.1)+DSTRS(NEL*2)+DSTRS<NEL*3) ) / 3. 

000342 

000 


DR ( 1 ) = SZZ 

000343 

000 


06(2) = 5RR 

0U0344 

000 


D£ 1 3) = STT 

000345 

000 


Ofi ( 4 ) = SZR 

000346 

000 


DFI = 0. 

000347 

000 


DFk =0. 

000348 

000 


00 200 I = 1*4 

000349 

000 


OF 1 = DFI + KL6U) * DSTKS ( NEL » T 1 

000350 

uou 

200 

CFk = OF* + OB ( I ) * CSTRS ( NEL #11 

000351 

000 


CFo = OElP * POP * XRS 

000352 

000 


ROW = -SmALLK / AlAMCA 

000353 

000 


DEN = (1 .+ETS/XMS) ** POW 

000354 

000 


DFk = DkN * ( DFK— o • *T J/P*PELP ) 

000355 

000 


OF = OFI - OFJ - bFK * ( 1 •- SR ALLK/ALAMDA ) 

000356 

000 


AS6 = XRS * P * (FOP-P) 

00U357 

000 


TTO = TO » 3. 

000358 

000 


S 108 A ( NEL ) = ASO 

000359 

000 


EP(NEL) = OF 

000360 

000 


WRITE <6.b2U) NEL * VO IOR * TT J » ASG* POP * OF . XMS » ETS 

000361 

000 

620 

FORMAT ( 1 5 * ' VOTOH=* *F10.4* » 3J=**F10.4** ASG= * * F10 . 4 * * P0=** 

00U362 

000 


IFIU.4*' 0F=«#E12.5*» XMS=* *E12.5* • ETS=**E12.5) 

000363 

000 

c 


00U364 

000 


IF(OF.LT.O. ) GO TO 764 

000365 

uou 

c 


000366 

000 


00 110 I = 1*4 

00U367 

000 


OB ( I ) - 0. 

000368 

uoo 


RO ( 1 ) = u. 

000369 

000 


00 110 J = 1*4 

000370 

Oou 


DR ( I ) = UR ( 1 ) + 0(1*0) * RLR(O) 

000371 

000 

110 

RD ( I ) = HO(I) + RLb(O) * 0(0*1) 

000372 

uoo 


OEh = 0. 

000373 

uoo 


DO 120 I = 1*4 

000374 

000 

120 

OEN = DEN + RLBU) * DBU) 

000375 

000 


OEi\i = OEN ♦ BR 

000376 

uoo 


DO 130 I = 1*4 

000377 

000 


DO 130 0=1*4 

000378 

000 

130 

DH.J) = 0(1*0) - UP(I) * RO(J> / DEN 

000379 

000 

764 

00 111 I = 1*4 

000380 

000 


00 111 0 = 1*4 

000381 

000 

111 

DRAT (NEL* I * J) = OU.J) 

000382 

oou 

600 

FORMAT ( 4E20 . 7 ) 

000383 

000 


RETURN. 

000384 

000 


ENO 

000385 

uoo 

lik.LT *SIH NASA*TPFS.STIFPl. * *113762121110 

000386 

000 


SUBROUTINE STIFFUNTAPE). 

0UU387 

000 

c---- 


OOU308 

00U389 

000 
o no 


PARAMETER NOOS=3UU*NELS=260*NF=200UO*MAX=600 

000390 

000 


COMMON /BLKO/ T IT lE ( 20 ) * INODE * NELER * NAPC * NBC * NINCR * IPRINT , EPSLON 

000391 

000 


COMMON /BLK1 / W(6) *H(6) *AR(4)*BR(4) *CR(4) * AZ ( 4 ) * BZ ( 4 ) * CZ ( 4 ) * 

000392 

oou 


* BN (4 ) * CN ( 4 ) *0N ( 4 ) . TYPEA ( 4 , 4 ) * TYPEB ( 4 » 4 ) * T YPEC ( 4 , 4 ) * TYPEE ( 4 * 4 ) * 

000393 

000 


* TYPEF (4*4) * T YFEG (4*4) r AO , BO *C0 *RT r RB *R A * RC * IC . JC * KC , LC * NEL 

000394 

oou 


COMMON /dLK6/ SIGH * SIGZ * S I GT * T AUZR *0(4*4) * STIFF (8*8) *KK( NELS* 8) * 

000395 

000 


1 R ( NODS ) * Z ( NODS ) *10T0IS<RAX) 

000396 

000 


COMMON ZBLK7 / DSThS ( NELS * 4 ) . ARM (NELS * 4 ) * AZR ( NELS * 4 ) *RTT (NELS) * 

000397 

000 


1 AOO(MELS) 


00004200 

00004300 

00004400 

00004500 

00004600 

00004700 

00004800 

00004900 

00005000 

00005100 

00005200 

00005300 

00005400 

00005500 

00005600 

00005700 

00005800 

00005900 

00006000 

00006100 

00006200 

00006300 

00006400 

00006500 

00006600 

00006700 

00006800 

00006900 

00007000 

00007100 

00007200 

00007300 

00007400 

00007500 

00007600 

00007700 

00007800 

00007900 

00008000 

00008100 

00008200 

00008300 

00008400 

00008500 

00000100 

■00000200 

00000300 

00000400 

00000500 

00000600 

00000700 

00000800 

00000900 

00001000 

00001100 

00001200 



1 


r- 


9 


000398 

000399 

000 

000 

c— 

COMMON /bLKA/M.N 


00001300 





000400 

000 

c 


* 

00001500 

000401 

000 

c 

DIFFERENT TYPE .OF INTEGRATIONS ARE 

PERFORMED 

AND SAVED ON FILE 00001600 

000402 

000 

c 

NO. 2 FOR LATER USE. 


00001700 

000403 

000 

c 



00001800 

000404 

000 


AZll) = Z<LC) - ZldC) 


00001900 

00040b 

000 


AZ ( 2 ) = Z(IC) - ZIKC) 


00002000 

00040ft 

000 


AZ13) = -AZ<1I 


00002100 

000407 

000 


AZ14)=-AZ12) 


00002200 

000406 

000 


GZ ( 1 ) = Z < KC ) - ZILC) 


00002300 

000409 

000 


eZ(2) = -BZ<1) 


00002400 

000410 

000 


ez<3) = ziua - zuc> 


00002500 

000411 

000 


BZ 14 ) = -BZ<3> 


00002600 

000412 

000 


CZll) = Z(KC) - ZldC) 


00002700 

000413 

000 


CZ 12 ) = ZUC) - Z ILC ) 


00002600 

000414 

000 


CZ13) = -CZ (2 ) 


00002900 

00041b 

000 


CZ14) = -C2(l) 


00003000 

000416 

000 


ARID = K ( dC ) - RILC) 


00003100 

000417 

000 


AR 1 2 ) = K(KC) - RUC) 


00003200 

000418 

000 


AR 1 3 ) = -ARID 


00003300 

000419 

000 


AR 1 4 ) = -AH (2) 


00003400 

000420 

000 


BR(1) = RlLC) - RIKC) 


00003500 

000421 

000 


BR (2 ) = -BR(1> 


00003600 

000422 

000 


BR 1 3 ) = RUC) - RldC) 


00003700 

000423 

000 


BR14) = -BR ( 3 ) 


00003600 

000424 . 

000 


CRll) = RldC) - RIKC) 


00003900 

000425 

000 


CR 12 ) = R(LC) - RUC) 


00004000 

000426 

000 


CR 13) = -CR 12) 


00004100 

000427 

000 


CR 14) = -CR(l) 


00004200 

000428 

000 


AO =-AR(3)*AZ(2) + ARl4)*AZll) 


00004300 

000429 

000 


BO — — BR l 2 ) *BZ l 4 ) + BR(3)*BZ<1) 


00004400 

000430 

000 


CO = CR (3) *CZ < 1 ) - CR ( 4 ) *CZ 12) 


00004500 

000431 

000 

c 



00004600 

000432 

000 


RT = RUC) + RldC) * RIKC) + R(LC) 


00004700 

000433 

000 


RA =-P(KC) + RUC) - RILC) + RldC) 


00004800 

000434 

000 


RB = RldC) - RUC) + RIKC) - RILC) 


00004900 

000435 

000 


RC =-RUC) + RldC) - R(KC) + R(LC) 


00005000 

000436 

000 


AOd(NEL) = AO 


00005100 

000437 

000 


RTT INEL) = RT 


00005200 

000438 

000 

c 



00005300 

000439 

000 


DO 200 ;M = 1.4 


00005400 

000440 

000 


ARMlNEL.M) = AKIM) 


00005500 

000441 

000 


AZMlNEL.M) = AZ(M) 


00005600 

0 0 04.^2 

000 


DO 200 N = 1.4 


00005700 

000443 

000 


CALL GAUSSU.AA) 


00005800 

000444 

oou 


CALL GAUSS 12 »BB) 


00005900 

000445 

000 


CALL GAUSS13.CC) 


00006000 

000446 

000 


CALL GAUSS14.EE) 


00006100 

000447 

000 


CALL GAUSS15.FF) 


00006200 

000448 

000 


CALL GAUSS (6 iGG ) 


00006300 

000449 

000 


TYREAIM.N) = AA 


00006400 

000450 

000 


TYPEBlM.N) = Bb 


00006500 

000451 

000 


TYPECIM.N) = CC. 


00006600 

000452 

000 


TYPEEIM.N) = EE 


00006700 

000453 

000 


TYPEFIM.N) = FF 


00006800 

000454 

000 


TYPEGIM.N) = GG 


00006900 


0UU455 
OOU456 
00U457 
00U456 
00U459 
00U46C 
0UU461 
00U462 
000463 
000464 
000466 
000466 
000467 
00U468 
000469 
000470 
000471 
000472 
000473 
000474 
000475 
000476 
000477 
000476 
000479 
000480 
000481 
000482 
00O483 
000484 
000485 
000486 
000487 
000488 
000489 
000490 
000491 
000492 
000493 
000494 
000495 
000496 
000497 
000498 
000499 
000500 
000501 
000502 
000503 
000504 
000505 
, 000506 
000507 
000508 
000509 
000510 
000511 


000 

OOO 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

000 

ooo 

000 

ooo 

ooo 

ooo 

ooo 

000 

ooo 

ooo 

000 

ooo 

ooo 

ooo 

Ooo 

ooo 

ooo 

ooo 

000 

ooo 

ooo 

Ooo 

ooo 

ooo 

ooo 

ooo 

ooo 

ooo 

ooo 

ooo 

ooo 


200 CONTINUE 00007000 

1*61 T£ (NT APE ) T Y PE A t TY PEB r T YPEC r TYPEE 1 T YPEF * TYPE6 00007100 

RETURN 00007200 

ENU 00007300 

i»ELT*SIh NASA*TPFS .STIFF2* « >113771121110 

SUBROUTINE STIFF2<nI*NTAPE) OOOOOlOO 

c 4 00000200 

PARAMETER NODS=30u * NELS=260 . NF=20000 *MAX=600 00000300 

C 00000400 

COMMON /6LK0/ T ITLE (20 ) * INODE * NELEM * N APC *NBC * NINCR r I PR I NT *EPSLON 00000500 
COMMON /bLK 1 / » <6> *H(6) * AR<4> *BR (4) *CR(4> , AZ(4) *BZ(4) *CZ<4) * 00000600 

* BN<4) *CN<4) .DN(4) * T YPEA (4 • 4 ) . T YPEB < 4 * 4 ) *TYPEC(4*4> *TYPEE(4*4) * 00000700 

* TTPEF (4*4) * T YPEG (4*4)* AO* 80 *CO*RT»RB*RA* RC * IC * JC • KC *LC * NEL 00000800 

COMMON /6LK6/ SIG«*SIG2*SIGT *TAU2R *0(4*4), STIFF (8*8>*KK( NELS* 8 ) * 00000900 

1 R (NCOS ) * Z ( NODS ) *TUTOIS(MAX) OPOOlOOO 

C 00001100 

c 00001200 

C FORM STIFFNESS MATRIX. 00001300 

c 00001400 

IF(Nl.NE.O) 00001500 

1READ (NTAPE) TYPEA*TYPEB*TYPEC* TYPEE*TYPEF,TYPEG 00001600 

. DO 200x1 = 1*4 00001700 

DO 200 J = 1*4 00001800 

c 00001900 

STIFF(I.J) = TYPtA(I*J)*0(l.l)/8.+TYPE8(l,J)*D(4,l)/ 8.+ 00002000 

1 TYPER (J*I)»Q(1*4)/ 8.+TYPEE(I*J)*D(4,4)/ 8. ' 00002100 

STIFF < J+4* 1 ) = TYPER( J* I ) *D(2, 1 ) /8.+TYPEC (J*I)«P(3*1) +TYPEA ( J* I ) * 00002200 

1 0(4*1)/ 8 «+T YPEE (U . I ) *0 (2 * 4 ) / 8.+T YPEF (<J*I> *0(3*4) *2+TYPEB < I * J) * 00002300 

2 D (4 *4 ) / 8. 00002400 

STIFF ( I »<J+4) = ST1FF(J+4,I) 00002500 

STIFF (1+4. J+4) = 1YPEE(I*vJ)*0(2.2)/8.+(TYPEB(I*J)+TYPEB<U*I>>* 00002600 

1 D( 2 , 4 ) / 8 .+TYPEA ( 1 * J ) *0(4*4) / 8 .+2 . * ( TYPEF (I * J ) +TYPEF ( J* I ) ) * 00002700 

2 0(2*3) + (TYPEC ( I * U) +TYPEC ( *)* I ) ) *0(3 *4 ) +TYPEG ( I , J) *D( 3* 3) 00002800 

200 CONTINUE 00002900 

RETURN 00003000 

ENU 00003100 

MELT *SIH NASA*TPFs.FRIC1N.-, *132675133010 

SUBROUTINE FRICTN ( 1C *KC*NEL,NFT , I SHEAR rVOIDI ) OOOOOlOO 

PARAMETER NOOS=3UO*NELS=26O,NF=2O0OO,MAX=6OO 00000200 

COMMON /BLK2/ ICMnODS , 2 > » I JKL (NELS * 4) , DEI *DE2 00000300 

COMMON /BLK4/ STRS (NELS *4 ) , DM AT < NELS* 4 * 4 ) , 0ELNM1 (MAX) . POP 00000400 

COMMON /BLK5/ DE (NtLS *4 *4 ) *SIGBA (NELS) rDSIGBA ( NELS) *DELL* YSTRS* 00000500 

1 FiNC*FN*ULOAO*FEt*PKAX*OLMAX 00000600 

COMMON- /BLK6/ SIGK,SI62%SIGT,TAUZR *D(4 ,4) , STIFF (8,8) ,KK (NELS, 8) , 00000700 

1 R (NODS ) • Z (NODS ) * T OTDI S ( MAX ) 00000800 

COMMON /6LK9/ PI » SMALL* *CK * BETA * PO *NFREE * NELST * IC ASE *NRIGd 00000900 

COMMON /BLK1U/ FRCK (20 ,8* 8 > * TR < 20 * 8 * 8 ) * XXL (20 ) *DZZ (20 > *DRR (20 > * 00001000 

1 POR (20 ) *DELTS *EC * XNUC *DELD OOOOllOO 

DIMENSION TS(8*8>- 00001200 

DZ = DZZ(NFT) 00001300 

DR = DRR.iNFT) 00001400 

PPI .4 3.14159 00001500 

RB = (R ( IC ) +R (KC ) ) A 2. 00001650 

PIR = PPI * RB / 3. 00001700 

BASE = R (KC ) - R 11C ) 00001800- 

HIGH. = Z(lt) - Z (KC) 00001900 

XL = (BASE+BASE + HTGH«HISH) **.5 00002000 



000512 

000 

XXl(nFT) = XL 

00002100 

0UU513 

000 

SINT = BASE / XL 

00002200 

0UU514 

000 

COST = HIGH / XL 

00002300 

000515 

000 

hRlTt(6ro00) NEL. XL. SINT. COST 

00002400 

0UU516 

000 

IF ( I SHEAR .NE . 0 ) 6g TO 160 

00002500 

000517 

000 

600 FORMAT!* NfcL .XL'. SINT. COST’. 15. 3F 12. 4) 

00002600 

00U518 

000 

CO 120 I = 1.8 

00002700 

000519 

000 

120 TR(NFT.I.I) = LOSI 

00002800 

000520 

000 

DO 130 I = 1.4 

00002900 

000521 

000 

J = 1 + 4 

00003000 

01)052? 

000 

TRtNF'T.I.J) = -SINT 

00003100 

001)523 

000 

130 TRtNF'T.J.I) = SIM 

00003200 

00U524 

000 

Dll = OZ 

00003300 

000525 

000 

Did = 0. , 

00003400 

000526 

000 

Odd = DR 

00003500 

000527 

000 

GO TO 17u 

00003600 

000528 

000 

160 CONTINUE 

00003700 

000529 

000 

NAmELlST/NAMEl/ SIGZZ.SIGRK.T AURZ* SIGN .SHEAR .TAUF.DI I »DZ 

00003800 

00O530 ' 

000 

II = ID(NFT.l) 

00003900 

000531 

000 

Jd = I0CNFT.2) 

00004000 

000532 

000 

SlG 22 = (STRSUI.l>+STRStJJ..l)> / 2. 

00004100 

000533 

000 

SIGRH = (STRSdl.^) +STKS l Jd . 2 ) ) / 2. 

00004200 

000534- 

000 

TAUR2 = tSTRSI 11 . 4) +STRSC Jd.4 ) ) / 2. 

00004300 

000535 

000 

SIGN = S1GZZ«C0ST*L0ST+S16RR*SINT*SINT-TaURZ*C0ST*SInT 

00004400 

000536 

000 

SHEAR = AHSISTKSlNtL.l) ) 

00004500 

000537 

000 

TAuF = CK ♦ S1GN»ulLTS 

00004600 

000536 

000 

Dll = DZ * (I.-SHLaP/TAUF) 

00004700 

000539 

000 

DZZ(NFT) = Dll 

00004800 

000540 

000 

Cld = 0. 

00004900 

000541 

000 

taR 1 Tt 1 6 . NAME 1 ) 

00005000 

000542 

000 

170 CONTINUE 

00005100 

000543 

000 

D I IF = Dll * P1R / XXL(NFT) 

00005200 

000544 

000 

DIdF = D1J * P1R / XXL ( NF T ) 

00005300 

000545 

000 

DJJF =, DdJ * PIP / XXL(NFT) 

00005400 

000546 

000 

DO 190 I = 1.4 

00005500 

000547 

000 

ST IFF < 1 . 1 ) = 2.+OUF 

00005600 

000548 

000 

STIFF (1+4. 1+4) = 2.*0JJF 

00005700 

000549 

000 

STIFF < I * 1+4 > = P.+UIJF 

00005800 

000550 

000 

190 ST IFF 1 1 . 9-1 ) = DUF 

00005900 

000551 

000 

CO 191 I = 1.2 

00006000 

000552 

000 

STIFF 1 1.1+2) r -DliF 

00006100 

000553 

000 

STIFF (1.5-1) = DliF 

00006200 

000554 

000 

STIFF U>-Z-l) = -2.*DIJF 

00006300 

000555 

ooo 

STIFF 1 1 . 1+.b ) = -OldF 

00006400 

000556 

000 

STIFF (T+4. 1+6) = -UJdF' 

00006500 

000557 

000 

ST IFF 1 1+4 . 9—1 ) = UJJF 

00006600 

000556 

000 

STIFF ( 1+2+ l+4> = -UIdF 

00006700 

000559 

ooo 

ST1FMI+2.9-I) - -2^*D1JF 

00006800 

000560 

ooo 

STIFF (2*1-1. 2*1) = -2.+DIIF 

00006900 

000561 

ooo 

191 ST IFF 12*1+3. 2*1+4 ) = -2.*DJJF : - 

00007000- 

000562 

ooo 

CO 210 I = 1.8- 

00007100 

.000563 

ooo 

DO 210 J = 1.8 

00007200 

000-564 

OOO 

210 ST IFF < J.-I) = STIFM1.JI . 

000073,00 

000565 

Ooo 

OMmT (NE t. 1.11 = Oil 

00007400 

000566 

ooo 

DRAT ( NEL » 1 . 2 ) = Did 

00007500 

000567 

ooo 

DMATU1EL.2.2) — Odd 

00007600 

000566 

ooo 

l*HITE.<6»630) NEL »D2 »GR *D1I .fllJ+Odd.XXKNFT ) 

00007700 


000569 

000 


630- 

•f ORMaT ( • NEL UZ OR oil DIJ DJJ LENGTH' . I5.6E12.5T 

00007800 

nuub7f.' 

000 



VsRITE (6.810) STIFF 

00007900 

0 0 u57 1 

000 


610 

FORMAT (8tl5.b) 

00008000 

OOU572 

000 



DO 140 I = 1.8 

00008100 

000573 

000 



DO 146 J = 1.8 

00008200 

00U57<- 

000 



SUfo = 0 • 

00008300 

(10057b 

000 



DO 14U K = 1.8 

00008400 

000576 

000 


144 

SON, = SUI* + TR(NFT.K.I) * STIFF(K.u) 

00008500 

00U577 

000 


140 

TS(I.J) = SUM 

00008600 

0O0S7P 

000 



DO 150 I = 1.8 

00008700 

00057s 

000 



CC 150 J = 1.8 

00008800 

000580 

000 



SUM =0. 

00008900 

000581 

000 



00 155 K = 1.8 

00009000 

000582 

000 


155 

SUM = SUM + TS(I.R) * TR(NFT.K.J) 

00009100 

000585 

000 


150 

ST IFF ( I. J) = SUM 

00009200 

000580 

000 



00 200 I = 1.8 

00009300 

00058.5 

000 



00 200 J = 1.8 

00009400 

00u58f? 

000 


20U 

FPCK(MFT.I.J) = ST IFF ( I »vl) 

00009500 

000587 

000 



NR LTE ( 6 » 610) STIFF 

00009600 

000588. 

ooo 



RETURN 

00009700 

000589 

000 



ENU 

00009800 

000590 

000 

L>ELT »SIH NA5A*TPFS.ASStMb». .114006121110 


QUUS91 

ooo 

C- 


SUbRCUTINE ASSEMblNEL.NET) 

ooonoioo 







00o593 

ooo 



PARAMETER NOOS=30U.NELS=260.NF=?un00»MAX=60n 

00000300 







000595 

ooo 



COMMON /bLKO/ T ITcE ( 20) . lNOHE .NFJLEM . NAPC . NBC . NINCR • IPRINT . EPSLON 

00000500 

000598 

ooo 



COMMON /bLN2/ lOtNOr.S.2) .IOKLINELS.4) .0E1.0E2 

00000600 

000597 

ooo 



COMMON /t)LK3/ XK(IMF) .APF(MAX) . IMAX » IBB . IHBI »LT. LAST 

00000700 

000598 

ooo 



COMMON /HLK6/ SIGR . S1G2 . SIGT . TAUZR » 0 (4 . 4 ) , ST IFF (8 . 8 > . KK ( NFLS. 8 ) . 

oooooBon 

000599 

ooo 



1 R (NoOS ) »Z (NODS ) .TOTDIS(MAX) 

00000900 

000600 

ooo 



COMMON /BLK1 0/ FRUK ( 20 . 8 . 8 ) . TR (20 .8 . 8 ) . XXL < 20 > .pZZ C 20 ) . DKR ( 20 ) . 

onooiooo 

000601 

ooo 

f . 


1 POR 120 ) . OELTS.EC . XNUC »OELD 

00001100 







000603 

ooo 

c 



00001300 

OOU604 

OOO 



NFREE = INODE * 2 

00001400 

000605 

ooo 



IFtNFT.LT.l.OR.NF'l .GT.NBC) GO TO 120 

00001500 

000606 

OOO 



DO 130 I = 1.8 

00001600 

000607 

ooo 



DO 130 J = 1.8 

00001700 

000608 

ooo 


130 

STIFFU.u) = FRCKINFT.1.J) 

00001800 

000609 

ooo 


120 

CONTINUE 

00001900 

000610 

000 

c 



00002000 

000611 

oou 



do no a = i .e 

00002100 

000612 

ooo 



II = KK(NEL.I) 

00002200 

000613 

ooo 



CO 110 J = 1.8 

00002300 

000614 

ooo 



Jd = KK(NEL.U) 

00002400 

000615 

ooo 



IFtll.EQ.O.OR.JJ.tO.O) GO TO 110 

00002500 

000616 

ooo 



IF(Il.LT.vJJ) GO TO 110 

00002600 

000617 

ooo 



IF(II.GT.IHBI) GO TO 104 

00002700 

000618 

ooo 



L = jJ + (11-1) *11/2 

00002800 

000619 

ooo 



GO TO 105 

0O002900 

000620 

ooo 


104 

L - «jJ 'M LT + (II-lHb) * IHQI 

00003000 

000621 

ooo 


105 

XK(U =-XK(L)-+ ST IFF ( I » J) 

00003100 

000622 

ooo 


110 

CONTINUE 

00003(?00 

000623 

ooo 



RETURN 

00003300 

000624 

ooo 



ENU 

00003400 

000625 

ooo 

MELT.SIH NA5A*TPF$ »D1SPL» ..114011121110 



000626 

nnn*27 

000 

noo 

C- 


subroutine dispumfree.ni. inode. incr> 

00000100 






-00000200 

000626 

000 



PARAMETER NODS=300 »NELS=260 ,NF=20000 »MAX=600 

00000300 

noo^2Q 

ooo 

c- 

■ d 

_ 







-00000400 

000630 

000 



COMMON /6LK2/ ID (NODS . 2 > . I JKL (NELS . 9 ) . DEl , DE2 

00000500 

000631 

000 



COMMON /6LK3/ XK (NF )» APF l MAX >» IMAX » IHB . IHBI . LT , LAST 

00000600 

000632 

000 



COMMON /BLK6/ S IGR » S IGZ , S IGT » TAUZR » D (9 . 9 ) .STIFF (6 .8 > .KK (NELS. 8) • 

00000700 

000633 

O0OA34 

000 

onn 

C- 


1 R(NODS) .Z(NODS) .TOTOIS(MAX) 

00000800 






-00000900 

000635 

000 

c 



.00001000 

000636 

000 



N — NFREE 

onooiloo 

00063? 

000 



ZTEST = 0.000001 

00001200 

000638 

000 

c 



00001300 

000639 

000 



DO 100 U = 1 .NFREt 

00001900 

000690 

000 



IF(J.GT.IHBI) 60 TO 106 

00001500 

000691 

000 



L = (J+l) * J / 2 

00001600 

000692 

000 



GO TO 109 

00001700 

000693 

000 


10 a 

L = LT + IHB * <J - IHBI) 

00001800 

000699 

000 


109 

XTEST = ABS(XKtL) 1 

00001900 

000695 

000 



IF(XTEST.LT. ZTEST) XK<L> =1. 

00002000 

000696 

000 


100 

•CONTINUE 

00002160 

000697 

000 

c 



00002200 

000698 

000 



DO 110 I = l.NFREE 

00002300 

000699 

000 


lio 

XK(LAST+I) = APF(l) 

00002900 

000650 

000 

c 



00002500 

000651 

000 

c 



00002600 

000652 

000 



CALL FACTOR (NFREt) 

00002700 

000653 

000 

c 



00002800 

000659 

000 

c 



00002900 

000655 

000 



CALL SOLTN(NFREE) 

00003000 

000656 

000 

c 



00003100 

000657 

000 

c 



00003200 

000658 

000 



WRITE (6 , 600 ) INCR.NI 

00003300 

000659 

000 



NN = 10 

00003900 

000660 

000 



DO 288 J = l.NN 

00003500 

000661 

000 



II = (J-l) * 2 + 1 + LAST 

00003600 

000662 

000 



JJ : II ♦ 1 

00003700 

000663 

000 


280 

WRITE (6,610) J»XK(1I) .XK(JJ>- 

00003800 

000669 

000 


600 

FORMAT!/// 20X, ’DISPLACEMENTS FOR CYCLE NO. ’ » I9//18X , *z - DISPL’. 

00003900 

000665 

000 



* 15X.’R - UISPL’.lbX.’LOAD INCREMENT STEP =»»I5//) 

00009000 

000666 

000 


610 

FORMAT (I7.2E20.7) 

00009100 

000667 

000 



RETURN 

00009200 

000668 

000 



END 

00009300 

000669 

000 

HELT»SIH NASA l *TPFS>. FACTOR. , ,119013121110 


000670 

000 



SUBROUTINE FACTOR (NFREt) 

00000100 

000671 

000 

c 


THIS SUBROUTINE PtHFQRMS FACTORING 

00000200 

000672 

000 



PARAMETER NODS=30U»NELS=260,NF=20000»MAX=600 

00000300 

000673 

000 



COMMON /BLK3/ XK(NF) .APflMAX), IMAX.IHB.IHBI.LT. LAST 

00000900 

000679 

000 



n = nfree 

00000500 

000675 

000 



IHbl = 1HB1 

00000680 

000676 

000 



DO 8 1=1, N 

00000700 

000677 

000 



IFII.GT.1HB1) GO 10 2 

00000800 

000678 

000 



K = 1 

00000900 

000679 

000 



M=K+(I-ll*l/2 

OOOOlOOO 

000680 

000 



GO TO 3 

00001100 

000681 

000 


2 

K=I-IH8l 

00003200 

000682 

000 



M=K-»LT+(1-IHW*IHB1 

00003300 



000683 

000 


3' 

•«J=I + 1HB1 

ononi4on 

000664 

000 



IFIJ.GT.N) GO TO 4 

00001500 

000685 

000 



JJ=I+IHB1 

00001600 

000686 

000 



GO TO 5 

00001700 

000687 

000 


4 

J J=N 

00001800 

000688 

000 


5 

6=0.0 

00001900 

000689 

000 



LA=1-1 

00002000 

000690 

000 



LB=I+1 

00002100 

000691 

000 



IF(LA.EQ.O) GO TO 6 

00002200 

000692 

000 



00 7 L=K » LA 

00002300 

000693 

000 



IF<L.GT.1HB1)G0 TO 50 

00002400 

000694. 

000 



J = < L+l ) *L/2 

00002500 

000695 

000 



GO TO 51 

00002600 

000696 

000 


50 

0 = LT + IHB*<L-lHol) 

00002700 

000697 

000 


51 

A = XK ( M ) 

00002800 

000698 

000 



R = B+A*A*XK(ol 

00002900 

000699 

000 


7 

M=M+1 

00003000 

000700 

000 


6 

A=XK<M) 

00003100 

000701 

000 



XK ( M ) = A-8 

00003200 

000702 

000 



|F ( I .EQ.N) GO TO a 

00003300 

000703 

000 



DO 9 J=Lt}> JJ 

00003400 

000704 

000 



SUM=0.0 

00003500 

000705 

000 



IF(J.GT.IHbl) GO TO 10 

00003600 

000706 

000 



K=1 

00003700 

000707 

000 



MM=K+< J-l)*J/2 

00003800 

000708 

000 



GO TO 11 

00003900 

000709 

000 


10 

K=J-1HB1 

00004000 

000710 

000 



KM=K+l.T+(J-IHB)*lhbl 

00004100 

000711 

000 


11 

IF(LA.EQ.O) GO TO 9 

00004200 

000712 

000 



1F< K.GT.LA) GO TO 9 

00004300 

000713 

000 



DO 12 JA=K*LA 

00004400 

000714 

000 



L=M-1+JA 

00004500 

000715 

000 



IF(JA.GT.IHBl) GO TO 13 

00004600 

000716 

000 



L1=(JA+1) * JA/2 

00004700 

000717 

000 



GO TO 14 

00004800 

000718 

000 


13 

L1=LT+IH8*< JA-lHBil 

00004900 

000719 

000 


14 

SUI»I=SIJM+XK(MM>*XML)*XK<L1> 

00005000 

000720 

000 


12 

WM=MM+1 

00005100 

000721 

000 


9 

XK(MM)=(XK (MM)“SUFt)/XKiM) 

00005200 

000722 

000 


8 

CONTINUE 

00005300 

000723 

000 



RETURN 

00005400 

000724 

000 



ENU 

00005500 

000725 

000 

4SELT rSIH NASA*TPFS.SOLTN« .» 114016121 1 10 


000726 

000 



SUbRGUT 1NE SOLTN ( nFREE fc 

00000100 

000727 

oou 



PARAMETER NOUS=3UU»NELS=260*NT=20000rMAX=600 

00000200 

000728 

000 



COMMON /8LK3/ XK (Nl ) » APF IMAX ) r IN AX . IFtB » IHBI » LT »L AST 

00000300 

000729 

000 

c 


THIS PORTION OF SubROUTlNE PERFORMS FORWARD-SUBSTITUTION 

00000400 

000730 

000 

c 



00000500 

000731 

ooo 



N = NFREt 

00000600 

000732 

000 



Ihbl. = I FIB 1- 

00000700 

000733 

000 



NF = LAST +1 

00000800 

000734 

000 

c 



00000900 

000735 

oou 


14 

DO I K = 2 »N 

00001000 

000736 

000 

c 



00001100 

000737 

000 



IFIK.GT.lHblJ GO Hi 2 

00001200 

000736 

000 



N = U 

00001.300 

000739 

ooo 



NN=K-1 • 

00001400 



000740 

000 



Ml=MN*K/2 

ooooison 

000741 

000 



GO TO 3 

00001600 

000742 

000 


2 

M=n-1HB 

00001700 

000743 

000 



MM— IHR1 

00001800 

000744 

000 



M 1=N* IHBl+LT 

0O001900 

000745 

000 


3 

SUM=0.0 

00002000 

000746 

000 



CO 4 L = 1 . MM 

00002100 

000747 

000 



LL=L+N 

00002200 

000748 

000 



JJ=LL+M1 

00002300 

000749 

000 



LL=LL*NF-1 

00002400 

000750 

000 


4 

SUM=SUM+XK ( UJ>*XK(LL) 

00002500 

000751 

000 


i 

XK(LL+1)=XK(LL+1 )-SUN 

00002600 

0007512 

ooo 



J = NF+N-1 

00002700 

000753 

000 

c 


THIS PORTION OF SUbRCUTlNE PERFORMS BACK-SURSTITUTION 

00002800 

000754 

000 



NF=NF+N-1 

00002900 

000755 

ooo 



XK (NF ) =XK (NF ) /XK (HAST ) 

00003000 

000756 

000 



DC 5 K=2»N 

00003100 

000757 

ooo 



L-N-K+l 

00003200 

000758 

ooo 



IF(L.GT.XHB1 ) GO TO 6 

00003300 

000759 

ooo 



I=L+(L-l)*L/2 

00003400 

000760 

ooo 



■ GO TO 7 

00003500 

000761 

ooo 


6 

I=L+(L-1HB)*IHB1+lT 

00003600 

000762 

ooo 


7 

IR=N-IHB 

00003700 

000763 

ooo 



IF(L.GT.IR) GO TO 8 

00003800 

000764 

ooo 



J=1HB1 

00003900 

000765 

ooo 



GO TO 9 

00004000 

000766 

000 


a 

J=K-1 

00004100 

000767 

ooo 


9 

SUM=0.0 

00004200 

000768 

ooo 



DO 10 Mzl iaJ 

00004300 

000769 

ooo 



MM=L+M 

00004400 

000770 

ooo 



IF(MM.GT.IHHI) GO TO U 

00004500 

000771 

ooo 



NN=L-MMM-l)*MM/2 

00004600 

0OU772 

ooo 



GO TO 12 

00004700 

000773 

ooo 


n 

NN=L+ (MM-IHR ) »lHbi+LT 

00004800 

000774 

ooo 


12 

MM=:NF-N+MM 

00004900 

000775 

ooo 


10 

SUM-SUM+XK (NN) *XK IMM ) 

00005000 

000776 

ooo 



MM— NF — N+L 

00005100 

000777 

ooo 


5 

XK(MM,)=XK(MM)/XK(i)-SUM 

00005200 

000778 

ooo 



RETURN 

00005300 

000779 

ooo 



END 

00005400 

000780 

ooo 

MELT »SIH NASA*TPF$. STRAIN. » .132705133010 


000781 

ooo 



SUBROUTINE STRAIN(NI» I NCR) 

00000100 

000782 


C- 









000703 

ooo 



PARAMETER N0US=30U » NELS=260 »NF=POOOO .MAX=600 

00000300 







000785 

ooo 



COMMON /BLKO/ TITLE ( 20 J . INODE .NELEM . NAPC .NBC > NINCR . IPRINT , EPSLON 

00000500 

000786 

ooo 



COMMON /BLK2/ ID ( NODS r 2 > « TUKL (NELS . 4 ) , DEI » DF2 

00000600 

000787 

ooo 



COMMON /BLK3/ XK (NF ) « APF (MAX ) .IMAX.IHP » iFtB I .LT .LAST 

00000700 

000788 

ooo 



COMMON /BLK4/ STH5 ( NELS. 4 ) » DM AT (NELS » 4 . 4 ) . DELNMK MAX ) . POP 

00000800 

000789 

ooo 



COMMON /bLK5/ DE ( Ntl.S .4 » 4 >-» SIGB A (NEL S ) » DS IGB A ( NELS ) .QELL'YSTKS. 

00000900 

000790 

ooo 



1 F INC » FN . ULOAD .FEL.PMAX. DLMAX 

nooniooo 

000791 ' 

ooo 



COMMON /BLK7/ USTKS ( NELS . 4 ) . ARM ( NELS > 4 > . A7M ( NELS r 4 > . RTT < NELS > . 

00001100 

000792 

ooo 



1 AOU(NELS) 

00001200 

000793 

ooo 



COMMON /BLK8/ 1NCNN . PDEPTH ( NELS ). VOTC I » ALAMDA .DEPTH . XMS 

00001.300 

000794 

ooo 



COMMON /6LK9/ PI » SMALLR# CK.RETA.PO.NEREE . NELST . I CASE rNR I6n 

00001400 

000795 

ooo 



COMMON /tJLKll/ VOID (NEL5J . DGAM ( NELS.4 > 

00001500 

000796 

ooo 



COMMON /BLK12/ SIGMX (NELS ). DEP (MELS ) »EP < NF.LS >« 0EGST2 (NELS ), ES 

00001600 


000797 

000 


- 

-CIXENSICN u<4) ,V<4) »EPS(4) ,SIG<4) 

00001700 

000796 

000 




00001800 

000799 

000 

L 



00001900 

000600 

000 



sura = o. 

0OQ02000 

000601 

000 



SU6B =■ o. 

00002100 

000802 

000 



IF l NfcC.l-.it . 0 ) CALL IN'TFaC (Del »Nl ) 

00002200 

000803 

000 



WP1T£<6,666> 

00002300 

000804 

000 


666 

F0KKAT(///10X, *STraINS A|\0 STRESSES FOR INCREMENT NO. =*»I4» 

00002400 

000805 

OQO 


* 

« • t-,0. OF CYCLES = •,14,//5X,*ELE6*»T14,*SIG-Z*.T30»*SI&-R*.T45. 

00002500 

000806 

OOU 


1 

t 'SIG-T* , T 60 , *TAU-aR* ,775, *AREA*»T9G» ’VOID RATIO* ) 

00002600 

000807 

oou 



CC 900 NlL = 1* NlLEM 

00002700 

000808 

000 



nft = NEL - nRIGD 

00002000 

000809 

0 00 



IF(NFT.G1 .O.ANO.NFT.LE.NbC) GO TO 900 

00 002900 

000810 

000 



oo loo i = i.4 

0O003000 

000811 

000 



IN = IJKL(NEL.I) 

00003100 

000812 

000 



II = <IN-1)*2 + L„ST + 1 

00003200 

000813 

000 



JJ = 11 + 1 

00003300 

000814 

000 



V<1) = XKUI) 

00003400 

000815 

000 


100 

U<1) = XK ( %JJ ) 

00003500 

000816 

000 

c 


FlUC STRAINS AnO SRtSSES AT THE CENTROID. 

00003600 

000817 

000 



E Z = 0. 

00003700 

000616 

000 



ER = 0. 

00003800 

000819 

oou 



Gfo = 0. 

00003900 

000820 

oou 



SUN. = 0. 

00004000 

000821 

000 



00 111 I r 1,4 

00004100 

000822 

000 



EZ = E2 ♦ ARM < NFL. 1) * VII) 

00004200 

000823 

000 



ER = ER + AZM(NEL.X) * Util 

00004300 

000824 

000 



66 = CM *■ ARM < NEL » 1 ) * U<I) + AZM(NFL.I) * V(I) 

00004400 

000825 

000 


111 

SUM = SUM ♦ UUI 

00004500 

000826 

000 

c 


COMPRESSION POSITIVE. 

00004600 

000827 

000 



EPS ( 1 ) = -EZ / AOo(NEL) 

00004700 

000628 

000 



EPS<2> = -ER / AOoU'EL) 

00004800 

0U0829 

000 



EPS ( 3 ) = -SUM / KIT(NEl) 

00004900 

00U830 

000 



EPS 14 ) = -GM / AO- (NEL) 

00005000 

000831 

000 

c 


EPS ( 1 ) = STRAIN IN 2 - DIRECTION. 

00005100 

000832 

000 

L 


EPS 1 2 ) = STRAIN IN R - DIRECTION. 

00005200 

000833 

000 

c 


EPS(3J = TANGENTIAL STRAIN. 

00005300 

000834 

000 

c 


EPS 1 4 ) = SHEAR STRAIN 

00005400 

000835 

oou 



CO 200 I s 1,4 

00005500 

000836 

000 



CGAMNEL.I) = EPSU) 

00005600 

000837 

000 



S IG ( T ) = 0. 

00005700 

000836 

000 



CO 200 J = 1,4 

00005300 

000839 

oou 



SIGtl) ,-= SIG(I) + DMAT (NEL » I , J) * tPS(J) 

00005900 

000840 

ono 


200 

CONTINUE 

00006000 

000841 

000 



IIi=NEL 

00006100 

00 08.42 

000 



DO 210 I = 1,4 

00006200 

000643 

oou 



DSTRS(II1.1)=S1G(1) 

00006300 

000844 

000 


210 

CONTINUE 

00006400 

000845 

000 



IF (NtL.LE. NRIGD) VjO TO 890 

00006500 

000846 

000 



CALL aREaA( IJKL(NtL.l) , IuKLlNFL ,2) , IJKL(NEL,3) , TJKL<NFL,4) , AKFA) 

00006600 

000847 

000 



RA)E = AREA / UELNM1 ( NEL ) 

00006700 

000846 

000 



VOID (NEL) = RATE * (l. + VOIDI) - 1. 

00006800 

000849 

000 


890 

CONTINUE 

00006900 

000850 

000 



ViRITE<6.600) NEL, S1G, AREA. VOID(NEL) 

00007000 

000851 

. 000 


900 

CONTINUE 

00007100 

000852 

000 

c 



00007200 

000853 

000 


600 

FORMAT (I9.7F15.6) 

00007300 


0008514 

000 

620 'FORMAT ( 110 * 2515.6) 

00007400 

000855 

000 

PE TURN 

00007500 

000856 

000 

END 

00007600 

000857 

000 

CiELT.SlH NASA*TPF$. INTP At, * , 13271 1133010 


000658 

000 

SUBROUTINE INTFAC(b'EI*NI) 

onoooioo 

000859 

000 

PArAMF TEH NOUS=30G*NELS=260*NF=20000*MAX=60n 

00000200 

000660 

000 

com ON /df-KC/ TITLt(20>*lNO0e*NELEM*MPC.N8e*NrNCR*NCYCL ,£PSLON 

00000300 

000661 

000 

COMMON. /6LK2/ ID ( NODS * 2 > * IUKL (PELS * 4 ) , OEl * DE2 

00000400 

000862 

000 

COMMON /bLK3/ XK(nF) .APF(MAX) » IMAX * IHR * IHBI *LT * L AST 

00000500 

000863 

000 

COMP, CM /tJLK4/ STKS (NELS . 4 ) * OP AT (NELS . 4 * 4 ) , OELNMHM AX > .POP 

00000600 

000868 

000 

COMMON /t)LN5/ DE(NELS.4.4) .SlGBA(NELS) * DSIGB A (NELS ) *DELL*YSTRS* 

00000700 

000865 

000 

1 F INC * FN * ULOAD *FEL*PMAX* DLMAX 

00000800 

000866 

000 

COMMON /BLK7/ UST«S(NELS*4) • ARP < NELS * 4 ) * A2M < NELS * 4 ) *RTT < NELS 1 * 

00000900 

000867 

000 

1 AOJ(NELS) 

oooniooo 

000868 

000 

C0MPUN/RLK6/ 1NCR.HDFFTh( NELS > * VOlDI *ALAMDA * DEPTH. PP 

onooi loo 

000869 

000 

COMPON /t>LK9/ HI*SPALLK.CK*BtTA*P0*NFREE*NELST.lCASE*NRI60 

0OU0120O 

000870 

000 

COMMON /6LK1U/ FRUK (20*8*8) *TR<20.a*8) » XXL (20) *DZZ<20) *0RO(20> * 

00001300 

000871 

000 

1 POR(?0> * DtLTS *EC • XNUC * DELO 

00001400 

000872 

000 

DIMENSION UG (8) *Ut- (8 ) 

00001500 

000673 

000 

DEI = 0. 

00001600 

000874 

000 

DO 900 NtL = 1*NBU 

00001700 

000875 

000 

III = NEL + NRiGu 

00001800 

000876 

000 

DO 100 I = 1»4 

00001900 

000877 

000 

II = IJKL(III.l) * 2 + LAST 

00002000 

000878 

000 

LGII) = XKUl-1) 

00002100 

000879 

000 

100 UGU+4) - XK(II) 

00002200 

000880 

000 

DO 110 1 = 1*8 

00002300 

000881 

000 

UL ( I ) = 0. 

00002400 

000882 

000 

DO 110 J = 1*8 

00002500 

000883 

000 

110 (JL ( I ) = UL ( I ) ♦ TK(NEL*I*J) » UG(J) 

00002600 

000884 

000 

WR ITt < 6 * 620 ) UG 

00002700 

000885 

000 

WRITE < 6 . 620 ) UL 

00002800 

000886 

000 

EZ - <-UL(l)*UC(2)+tlL<3)-UL<4)) / <2.*XXL (NEL > ) 

00002900 

000887 

000 

ER = (~UL<5)+UL(6)+UL<7)-UL<8>) / (2.»XXL (NEL) > 

00003000 

000888 

000 

SIGZ = Dm AT (111*1*1) * EZ + DMAT (111*1*2) * ER 

00003100 

000889 

000 

SIGR : DMAT (111*1*2) * EZ + DMAT ( 1 1 1 * 2 * 2 ) * ER 

00003200 

000890 

000 

SIGZ = -SIGZ 

00003300 

000891 

000 

CSTRS(Hi.l) = SIGZ 

00003400 

000892 

000 

DSTRS ( 1 1 1 *,2 ) = SIGN 

00003500 

000893 

000 

NAMELIST /NAME2/ I i 1 *SIG?*SlGR 

00003600 

000894 

000 

WPlTfc(6.NAME2) 

00003700 

000895 

000 

900 CONTINUE 

0O003800 

000896 

Onu 

620 FORMAT (HF 15.5) 

00003900 

000897 

000 

PE (URN • 

00004000 

000898 

000 

ENU • 

00004100 

000899 

000 

l»ELT*SIH NASA*TPFS. SETUP* * *114055121 110 


000900 

000 

SUbROUTINE SETUP 

oooooloo 

000901 

000 

COMMON /6LK1/ W < 6 ) *H ( 6 ) * AR ( 4 ) *BR ( 4 ) * CR ( 4 ) , AZ ( 4 ) , BZ < 4 ) *CZ ( U ) . 

00000200 

000902 

000 

* BN ( 4 ) * Cn ( 4 ) *Dn(4) *TYPEA(4*4)*TTPEB(4,4) * TVPEC ( 4 * 4 ) » T YPEE ( 4 * 4 ) * 

00000300 

000903 

000 

* TYPEF<4*4)*TYREG(4,4)*A0.B0*C0*RT*RB*RA*RC*IC.JC.KC,LC.NEL 

00000400 

000904 

000 

W(l) - .1713244924 

00000500 

000905 

000 

W ( 2 > = .3607615730 

00000600 

000906 

000 

w ( 3 ) = .4679139346 

00000.700 

000907 

000 

W(4) = W(3) 

OOOOOLOO 

000908 

000 

, W ( 5 ) = W < 2 > 

00000900 

000909 

000 

W<6). = W(l) 

oooniooo 

000910 

000 

Hill = .9324695142 

00001100 


000911 

000 

M2) = .6612093665 

00001200 

000912 

000 

M3) = .2386191661 

0000)300 

000913 

000 

M 4 ) = — H ( 3 ) 

00001400 

000914 

000 

H(5) = -H<2) 

00001500 

000915 

000 

M 6 ) — — H ( 1 ) 

00001600 

000916 

000 

BN ( 1 ) = -1. 

00001700 

000917 

000 

BN ( 2 ) = 1. 

00001800 

000918 

000 

BN ( 3 ) = 1. 

00001900 

000919 

000 

BN ( 4 ) = -1. 

00002000 

000920 

000 

CN ( 1 ) = 1. 

00002100 

000921' 

000 

CN ( 2 ) = 1. 

00002200 

000922 

000 

CN ( 3 ) = -1. 

00002300 

000923 

ooo 

CN (4 ) = -1. 

00002400 

000924 

000 

ON ( 1 ) = -1. 

00002500 

000925 

000 

ON ( 2 ) = 1. 

00002600 

000926 

000 

ON { 3 ) = -1. 

00002700 

000927 

ooo 

ON (4 ) = 1. 

00002800 

000928 

000 

RETURN 

00002900 

000929 

ooo 

END 

00003000 

000930 

ooo 

liELT »SIH NASA*TPFS.ELASTL. » *114060121110 

** 

000931 

ooo 

SUBROUTINE ELASTC (U • NR I GO » £ » XNU ) 

oooooion 

000932 

ooo 

DIMENSION 0 ( 4 #4 ) 

00000200 

000933 

ooo 

WRITE (6.600 ) NRI6U.E « XNU 

00000300 

000934 

000 

CONST = E*XNU / { U.+XNU)*(1.-XNU*2.) > 

00000400 

000935 

ooo 

SHEAR = E / (2.»(1.+XNU) ) 

00000500 

000936 

ooo 

0(1.1) = CONST + SHEAR*2. 

00000600 

000937 

ooo 

0(2.2) = 0(1.1) 

00000700 

000938 

ooo 

D (3.3 ) = 0(1.1) 

00000800 

000939 

ooo 

D ( 4 . 4 ) r SHEAR 

00000900 

000940 

ooo 

0(1.2) = CONST 

00001000 

000941 

ooo 

0(1.3) = CONST 

00001100 

000942 

ooo 

0(2.3) = CONST 

00001200 

000943 

ooo 

DO luO I = 1.4 

00001300 

000944 

ooo 

CO 100 J = 1.4 

00001400 

000945 

ooo 

100 D(J.I) = 0(1. U) 

00001500 

000946 

Ooo 

600 FORMAT ( // ' FOR FIRST* .14,* ELEMENTS. THE FOLLOWING MATERIAL PRO 

00001600 

000947 

ooo 

* PERTIES ARE USED TO FORM ELASTIC MATRIX (0)'/* E ='.F20,7. 

00001700 

000948 

OOO 

* • XNU -* »F1U. 3//) 

00001800 

000949 

OOO 

RETURN 

00001900 

000950 

000 

END 

00002000 

000951 

Ooo 

ftELT.SIH NASA*TPFS. GAUSS... 114063121 110 


000952 

OOO 

SUBROUTINE GAUSS UT.AA) 

00000100 

000953 

OOO 

COMMON /dLKl/ W (6 ) »H (6) . AR (4) .BR (4 ) .CR ( 4 ) . AZ (4 ) , B? (4 ) »CZ ( 4 ) . 

00000200 

000954 

OOO 

* BN (4 ) »CN(4> .ON (4) .TTPEA (4.4) .TVPEB (4.4) .TTPEC (4 ,4 ) ► T YPEE (4.4) . 

00000300 

000955 

ooo 

* TYPEF(4.4) .TYPEG(4.4) .AO.BO.CO.RT.RB.RA.RC.IC.JC.KC.LC.NFL 

00000400 

000956 

ooo 

TWOPI = 6.28318531 

00000500 

000957 

ooo 

IPT = 6 

00000600 

000958 

ooo 

AA = 0. 

00000700 

000959 

ooo 

00 100 I = 1 .IPT 

00000800 

000960 

Ooo 

X '= Ml) 

00000900 

000961 

.000 

DO 100 ,J = '1 * IPT 

onooioon 

000962 

ooo 

Y = H(J) 

00001100 

000963 

ooo 

AA = AA * w(l) * *(J) * F(X.T.IT) 

00001200 

000964 

ooo 

100 CONTINUE 

00001300 

000965 

ooo 

AA = AA * TWOPI 

00001400 

000966 

ooo 

RETURN 

00001500 

000967 

000 

END 

00001600 


000968 

000 

MELT » S IH NASA* I PF5.F t t > 114066121110 



000969 

000 


FUNCTION F(X»Y»IT> 


00000100 

000970 

000 


COMMON /bLKl/ A (6) fH(6) . AR(4) »PR (4) .CP (4 ) , A> ( 4 ) ,HZ(4 ) >CZ<U ) , 


00000200 

000971 

000 


* RN<4) >CN(4) »DN<4) *TYPEA(4#4) *TYPEb<4.4) .TYPtC(4*4) »TYPEE(4> 

4) , 

00000300 

000972 

ooo 


* T YPEF ( 4 > 4 ) »TYFEG(4.4) » A(j # t>0 » CO » RT » Rb » RA » RC r IC « JC » KC » LC » NFL 


00000400 

000973 

000 


COMMON /bLRA/M»N 


00000500 

000974 

000 

c 

X STANDS FOR ZhAl In /HAI - EITA COORD. 


00000600 

000975 

ooo 

c 

Y STANDS FOR tlT A In ZHAI - EITA COORD. 


00000700 

000976 

ooo 

c 

Fc = UET . OF JACOB I • 


00000800 

000977 

000 

c 

FC = IN) * (R) 


00000900 

000978 

ooo 


FB = AO + 60 * X + CO * 1 


onooiooo 

000979 

OGO 


FC = < RT + RB * X + RA * Y + RC * X * Y) / 4. 


00001100 

000980 

ooo 


60 TO (10»20«30»4U«50#60) * IT 


00001200 

000981 

ooo 

10 

F = (AR(M)+BR(M>*a+CR(M>*Y) * (AR(N>+PR(N>*X+CR(N)*Y> /Fb 


00001300 

000982 

ooo 


F = F * FC 


00001400 

000983 

OOO 


RETURN 


00001500 

000984 

oeo 

20 

F = (AZ(M)+RZ(M)*X+CZ<M)*Y) * (AR(N)+PR(N)*X + CR(N)*Y.) /Fb 


00001600 

000985 

ooo 


F = F * FC 


00001700 

000986 

ooo 


RETURN 


00001800 

000987 

ooo 

30 

F = (l.+dN(M)*X+Ci\,(v)*Y+ON(M)*X*Y) * ( AR ( N > +bR (N ) *X+CR ( N > * Y ) 

/ 32 

.00001900 

000988 

000 


RETURN 


00002000 

000989 

ooo 

40 

F = (AZ(M>+PZ(M>*X+CZ(M>*Y> * t AZ ( N ) +PZ < N ) *X+CZ ( N > *Y ) /Fb 


00002100 

000990 

000 


F = F * FC 


00002200 

000991 

ooo 


RETURN 


00002300 

000992 

ooo 

50 

F = ( 1 • +bN ( M ) *X+CN ( N ) *Y+UN ( M ) *V*Y > * ( AZ (N ) +BZ t N ) *X+CZ ( N ) * Y ) 

/ 64 

.00002400 

000993 

OOO 


RETURN 


00002500 

000994 

ooo 

60 

F = (l.+dN(M)*X+CN(M)*Y+DN(M)*X*Y) * 


0n002600 

000995 

ooo 


1 <1.+RN(n)*X+CN<N)*Y+0N(N)*X*Y> * FP / <128. *FC) 


00002700 

000996 

OOO 


RETURN 


00002800 

000997 

ooo 


END 


00002900 

000998 

ooo 

WELT • SIH NaSA*TPF1>. AKEAA» • » 114U72121 1 1 0 



000999 

ooo 


SUBROUTINE AKEAA (1C,JC»KC.LC. AREA) 


onoooloo 

001000 

ooo 


PARAMETER NOUS=3uu.NELS=2bO.NF=20000»MAX=600 


00000200 

001001 

ooo 


COMMON /dLKb/ SIGK.SIGZ.SIGT.TA(IZR.D(4.4) ,STIFF(8.8> »KK(NELS»8) . 

00000300 

001002 

ooo 


1 R ( NCOS ) .7 (NODS) . lOTUIS(MAX) 


00000400 

001003 

ooo 


A I = (R< JC)-R< IC) ) * (Z(LC)-Z(IC) ) - (R(LC)-R(IC) ) * <Z<JC)- 

Z ( IC) >00000500 

001004 

ooo 


AJ = (R(KC)-R(JC) ) * (Z<LC)-Z< JC) ) - (R(LC)-R(JC) ) * (Z<KC)- 

Z<JC) >00000600 

001005 

ooo 


IF(Al.LT.O) AI = -Al 


00000700 

001006 

ooo 


IF(AJ.LT.O) AJ = -AJ 


00000800 

001007 

ooo 


AREA = (AT + AJ) / 2. 


00000900 

noiooa 

ooo 


RETURN 


onomoon 

001009 

ooo 


ENU V 


00001100 

001010 

ooo 

loiELT *SIH NASA*TPFl.PTLCMUn» » 114074121 1 10 



001011 

ooo 


SUBROUTINE PTLOAQ (NAPC » ULOAD ) 


ononoioo 

001012 

ooo 


PARAMETER NODS=3uU»NELS=260.NF=200UO»MAX=600 


00000200 

001013 

ooo 


COMMON /bLK3/ XK (NF ). APF (MAX ) . IMAX r IHP » IHBI » LT . LAST 


00000300 

001014 

ooo 

C 



00000400 

001015 

ooo 

c 



00000500 

001016 

ooo 

c 

get concentrated load if there is any 


00000600 

001017 

OOO 


WRITE (6»b77 ) 


00000700 

001018 

ooo 


DO 920 NC = l.NAPC 


00000800 

001019 

OOO 


REAU (5*540) NODE » PZ • PR 


00000900 

001020 

ooo 


WRITF(6»688) MODE » PZ » PR 


ononioon 

001021 

ooo 


II = (NOUE-1) * 2 


00001100 

001022 

ooo 


APF<1I+1) = PZ 


00001200 

001023 

ooo 


APF (11+2) = PR 


00001300 

001024 

ooo 


ULOAD = PZ 


00001400 


01)1025 

000 


920 CONTINUE 

ononisoo 

00102a 

000 


54 U FOKNhT ( I s #2Fl 5* 5 ) 

00001600 

001027 

000 


677 FOhM.AT (///10X, * APPLIED C. LOAD • //5X ,• NODE *. I OX .* FORCE TO 7*,5X»* 

00001700 

001020 

000 


* force to k*//) 

00001800 

001020 

000 


66,6 FORMAT (5X,i5>2t5X,L12. 4) ) 

00001900 

0U1030 

000 


RE 1 URN 

00002000 

001031 

000 


ENU 

00002100 

001032 

oou 

LELT »SIH NASA*TR6'S.F:0LU„U, , ,114076121110 


001033 

oou 


SUBROUTINE EuLUAo iULOAU) 

ononnioo 

001U30 

ooo 

c 


00000200 

001035 

000 

c 

THIS SUoROliTlNc. TS TO CALCULATE THF EQUIVALENT NODAL LOADS. 

00000300 

001036 

000 

L 


00000400 

001037 

oou 


PARAMETER NOUS=30G , NELS=26U , NE=?0000 , 4>AX=600 

00000500 

001036 

oou 


COM.MON /bl.K?/ 1DIWUI'S,2> »IoKHNELS» 4> ,0Fl,l)E2 

00000600 

001039 

000 


COMMON /bLK3/ XK (n 6 ) ,APF (MAY) , TMAX, IhB, IHBI ,LT, LAST 

00000700 

OU1U40 

oou 


COMMON /BLK6/ S I UK » S 1GZ > SIGT » TAUZR , D 1 4 , 4 ) ,STIFF(6,8) »KK1NFLS»8> , 

00000800 

001041 

ooo 


1 RINCDS) ,Z(|10DS) , fOTDIS(MAX) 

00000900 

001042 

000 


PI = 3.1415926 

oooomoo 

001043 

oou 


XXI = 0. 

onooiloo 

001044 

ooo 


XXU : 0. 

00001200 

001045 

OOO 


RE AO 15 * 500 ) NLDF.L , ULO AU 

00001300 

001046 

ooo 


WRITt(6,600) NLDti. 

oooniuon 

001047 

ooo 


CO 2uO 1 = l.NLDtL 

00001500 

001048 

ooo 


READ (5, 510} N0L,NL6T#NRHT 

00001600 

001049 

ooo 


KC = NRHT 

ononi7on 

001050 

ooo 


LC = NLFT 

00001800 

001051 

ooo 


EGL = (RlKC)**2-RlLC)**2)*PI*UL0AD/2. 

00001900 

001052 

ooo 


11 = (KC-1) *2+1 

00002000 

001053 

oou 


JJ = (LC-1) *2+1 

00002100 

001054 

ooo 


APR ( 1 I ) = APF(XI) + E0L 

00002200 

001055 

ooo 


APF(viJ) = aPF(UJ) + EQL 

00002300 

001056 

ooo 


200 V.RITt(6,b20) NOI. , uC , KC , FGL fULOAD 

00002400 

001057 

ooo 


620 F0HMAT(I6,5X,214,262U.7) 

00002500 

001056 

ooo 


500 FORMAT ( 15,620.9) 

00002600 

001059 

ooo 


510 FORMAT ( 315 ) 

00002700 

001060 

UlJO 


600 format i /// * calculation of equivalent nodal load*/ 

00002800 

001061 

MOO 


* • NUMBER uF LOADED ELEMFNTS = ',14/ 

00002900 

001062 

ooo 


* ' ELEM NO LOADED NODE EQUIV LOAD GIVEN U.LOAO*/) 

00003000 

001063 

ooo 


610 FOKMAT12ilO,2F15.7 ) 

00003100 

001064 

ooo 


RETURN ' 

00003200 

001065 

ooo 


ENU 

00003300 

001066 

Oou 

Wfc.LT »SIH NASA* TPHS.ZtPO, ,,114100121110 


001067 

uoo 


SUBROUTINE ZERO 1 A , N , N ) 

00000100 

001066 

oou 


dimension aid 

00000200 

001069 

ooo 


K = N * Mi 

00000300 

001070 

oou 


do ion i = i,k 

00000400 

001071 

ooo 


100 All) = 0. 

00000500 

001072 

Ooo 


RETURN 

00000600 

001073 

ooo 


END 

00000700 

001074 

ooo 

l»M AP , 1 X 


001075 

uoo 

tXPT 


001076 

ooo 


CONE PtNETROME TER 


001077 

ooo 


90 71 3 1UU 52 


0011)76 

ooo 


25. 


001079 

uoo 




001080 

oou 


7000000. 0.3 10. 0.45 


noioai 

000 
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APPENDIX 2 
DATA INPUT FORMAT 

Card 1: FORMAT (20A4) 

TITLE - Title of the problem 
Card 2: FORMAT (1015) 

(1) INODE - No. of nodes 

(2) NELEM - No. of elements 

(3) NAPC - No. of applied point load 

(4) NBC - No. of interface element 

(5) NINCR - No. of load increment 

(6) NCYCL - a dummy 

(7) ICASE = 0 for plastic analysis 

(8) NRIGD - No. of rigid element 

(9) NULOAD - No. of uniformly loaded element 

Card 3: FORMAT (8F10.4) 

(1) YSTRS - a dummy 

(2) DELL 

(3) ZETA - " 

(4) PMAX - Maximum load one wants to apply 

(5) DLMAX - a dummy 

Card 4: FORMAT (4F20.5) 

t (1) DZI - Shear modulus for interface element 

(2) DRI - Rotational modulus for interface 
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Card 


Card 


Card 


Card 


5: 

FORMAT 

(4F20.5) 

(1) 

EC 

- Modulus of elasticity for rigid element 1 

(2) 

XNUC 

- Poisson's ratio for rigid element 

(3) 

ES 

- Modulus of elasticity for soil 

(4) 

XNUS 

- Poisson's ratio for soil 

6: 

FORMAT 

(8F10.4) 

(1) 

PI 

- Angle of friction in degree 

(2) 

SMALLK 

- Swelling index 

(3) 

XI 

- Adhesion (for interface element) 

(4) 

VOIDI 

- Initial void ratio 

(5) 

PO 

- Initial density 

(6) 

DEPTH 

- Maximum depth of soil 

(7) 

ALAMDA 

- Compression index 

(8) 

EPSLON 

- Angle of friction for interface element 

7 = 

FORMAT 

(2F10.4, 215) 

(1) 

Z(I) 

- Z - Coordinate value (downward positive) 

(2) 

R(I) 

- R - Coordinate value 

(3) 

IZ 

= 0 if free to Z-direction 



= 1 if note 

(4) 

IR 

= 0 if free to R-direction 


= 1 if note 

Repeat INODE times in the order of node number 
8: FORMAT (415) 

4 node numbers of an element in counter-clockwise. Repeat 
NELEM times in the order of element number. Ordering of 


element should be: 
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(1) Rigid element (2) interface element 
(3) soil element 
*Ca rd ( s ) 9: FORMAT (1015) 

(1) ID(T.,1) - Element number left to I— interface element 

(2) ID(I,2) - Element number right to I— interface element 

Repeat NBC/5 times with 5 sets of data on one card 

* 

Not required if NBC=0 
Card (s) 10: FORMAT (15, 2F15.6) 

(1) NODE - Node number with point load 

(2) PZ - Z-component 

(3) PR - R-component 

Repeat NAPC times 

Not required if NAPC=0 
Card 11: FORMAT (15, F20.9) 

(1) NLDEL - No. of uniformly loaded elements 

(2) ULOAD - Load intensity (compression is positive) 

Card (s) 12: FORMAT (315) 

(1) NOL - Loaded element number 

(2) NLFT - Node No. at left 

(3) NRHT - Node No. at right 
Repeat NLDEL times 

** Not reguired if NULOADsO 


r 


APPENDIX 3 


FLOW CHART 



A 


















Form stiffness matrix 
for rigid element 


8 
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APPENDIX 4 

SUBROUTINE ORGANIZATION CHART 
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Subroutine 
Name 

SETUP 

INPUT 

ZERO 

DMATRX 

STIFF 1 

GAUSS 

F 

FRICTN 

AREAA 

ELASTC 

ASSEMB 

PTLOAD 

EQLOAD 

DISPL 

FACTOR 

SOLTN 


APPENDIX 5 

DESCRIPTIONS OF SUBROUTINES 


DESCRIPTIONS 


Assigns necessary constants for integrations 

Reads and writes input data, and gets ready for the 
first linear elastic solution. 

Clear a given matrix. 

Computes D^ ^ and forms D = D( ^ + D( j 
if dF> 0 Tor a given elementT 

Forms submatrices of stiffness matrices 

Integrates by Gaussian quadrature 

Gives functions to be integrated 

Updates interface moduli and forms interface element 
stiffness . 

Computes cross sectional area of an element. 

Forms elastic matrix D( 0 ) 

Assembles stiffness matrices into global stiffness 
matrix applying boundary conditions'. 

Reads and writes applied point load if there is any. 

Reads and writes applied uniform load and computes 
equivalent nodal forces. 

Calls FACTOR and SOLTN, and writes du for first 10 nodes. 

Factors the given simultaneous eqns. in one dimensional 
’ ’ array. 

Backward substituion is performed to give a set of 
1 solutions. 


STRAIN 


Computes incremental strains and stresses. Also 
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INTFAC 

STIFF2 


i 


computes new void ratio. 

Computes incremental stresses for interface element 
if there is any. 

Forms stiffness matrix using D. 
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PART II 

DYNAMICS OF WHEEL- SOIL INTERACTION 
II -1. INTRODUCTION 

Deformations and stresses of soil media under a moving wheel 
are complex phenomena. A rate-dependent inelastic behavior associated 
with inertia effects must be considered. Somewhat simplified analyses 
have been reported by various authors. Earlier contributions to the 
wheel-soil interaction by Bekker were followed by Micklethwaite £2], 
Evans [3], Uffelmann [4], and Bekker [5]. Rigorous experimental and 
theoretical studies on this subject have also been reported by 
Onaffeko and Reece [6], Wong and Reece [7,8]. Yong and Webb [9] and 
Schuring [10] studied energy dissipation in soil-wheel interaction 
from the viewpoint of viscoplasticity. Windisch and Yong [11] further 
examined the strain-rate phenomena and presented a method of computing 
soil displacements and strain rates from the experiment-based "marker 
position". In contrast to these studies, Perumpral, Liljedahl and 
Perloff [12] used the finite element method to calculate stresses and 
deformations due to a rigid wheel interaction. They used variable 
modulus of elasticity determined from the stress-strain curve of the 
triaxial tests but ignored the effects of inertia and rate-dependency. 

Elsamny and Ghobarah [13] studied the stress field in the soil 
mass under the loading of a rigid cylindrical wheel on the verge of 


spinning. However, the fact that the kinematic characteristics of the 
wheel and the velocity boundary conditions on the wheel-soil interface 
is ignored has been criticized by Wong [14] . More recently, Kloc [15] 
presented analytical formulations on mechanical interaction of a driven 
roller on soil slopes. In this study, a gravitating cohesive-frictional 
soil was considered with Kotter's quasi-static equilibrium equations 
applied to a plastic stress configuration (Mohr-Coulomb criteria) sat- 
isfying Shield's velocity conditions along the characteristic lines. 
Energy dissipation was not considered in this study. 

In the present study we propose a rational approach in which 
the rate-dependent inelastic properties together with effects of inertia 
are adequately taken into account. Equilibrium conditions for wheel- 
soil interaction reported by Onaffeko and Reece [6] and Wong and Reece 
[7] are used to obtain radial and tangential stresses at the interface. 
Although the nonisothermal conditions may be considered without special 
difficulties in the framework of continuum mechanics and irreversible 
thermodynamic process, the present study is limited to an isothermal 
condition. The Mohr-Coulomb failure criterion appears to dominate most 
of the wheel-soil interaction studies. However, in view of the fact that 
the soil behaves as a strain-hardening material, in general, rather than 
perfectly plastic or rigid plastic material, we will overcome such defi- 
ciency by using the concept of critical state soil mechanics. 

In what follows we make use of the internal state variable 
approaches of Coleman and Gurtin [16] and Perzyna and Wojno [17] 



11-54 


However, a basic difference from their approach is introduced in the 
present study such that the free energy functional containing inelastic 
behavior is not considered smooth for its entire domain of histories. 
Rather, we assume a form of discretized free energy as a function of 
elastic strains, plastic strains and internal or hidden variables of 
incremental quantity considered to be valid only for a small time 
interval or a fraction of loading increments. Here the hidden varia- 
bles may represent a viscous or physico chemital' behavior, properties 
other than what is commonly known as "elastic" and "plastic". Once 
the form of incremental free energy containing all nonlinear functions 
is prescribed for a small time interval, then the superposition of 
these nonlinear terms is permissible. Namely, the plastic material 
kernel may be calculated from the independent viscoelastic responses 
within this small time interval. Thus the histories can be carried 
over from one time increment to another until desired histories are com- 
pleted. This will be accomplished by a suitable difference operator. 

To represent inelastic behavior of soil we use the concept of 
critical state [18j and yield surface of Roscoe and Burland [1? ]. 

A derivation of the plastic tangent stiffness matrix based on this 
theory in the context of incremental theory of plasticity and its fi- 
nite element applications were presented in Part I of this report. It 
should be noted that the particular internal state variable approach 
used here in conjunction with incremental free energy expression leads 
to a valid coupling of the completely independent plasticity theory and 
the rate dependent hidden variables. 

Numerical examples are presented to demonstrate effectiveness of 
the present method. The well-known finite element method (t.O,£l j is 
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utilized in the computation. 

II - 2. BALANCE OF ENERGY AND LINEAR MOMENTUM 

We record here the principle of conservation of energy which 
states that the time rate of change of the kinetic energy k plus the 
internal energy U is equal to R, the mechanical power on the system. 


£ + U = R 

Here the superposed dot indicates a time rate, and 

1 


k = ^ I oVjVjdV 

*V 


I: 

I 

R*| pF^VjdW- I s^VjDjdA 

* v «/a 


u =| p p dv 

V 


( 1 ) 

( 2 ) 

(3) 

(4a) 


in which p is the density, v, is the velocity component; G is the 

. i i J 

internal energy density; F is the body force; s is the surface 

traction; and n t is the unit normal to the surface. Using the Green - 

Gauss theorem, (4a) becomes 


J 


R = I (pF* Vj + + a> j v j )dV 


(4b) 


Now, inserting (2) and (4b) into (1) yields 


Ji, 


[(a M + pF - pa )Vj - p G + <j V j>i3^ v = ® 


(5) 


For the principle of balance of linear momentum to hold and for 
arbitrary volumes we must have 


‘J L i - A 

C7> J + pF " P a - “0 


( 6 ) 
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and 


. i J _ 1J . 

pc ~ a v j m “ a V t j 


(7) 


i i 


Here and Y t ; are the stress tensor and strain tensor; the comma 
denotes ordinary differentiation; and a i is the acceleration. It 
should be noted that equations (2) through (7) refer to rectangular 
cartesian coordinates. We regard (7) as the balance of energy. 


II - 3. INCREMENTAL FREE ENERGY FUNCTIONS 

In view of the earlier discussion our objective is to propose 

a form of free energy functions in incremental quantity such that 

the non-smooth or inelastic strains may be included for a small time 

interval At. For isothermal conditions, the incremental free energy 

1 J 

$ (At) and stresses CT (At) are assumed to be functions of Incremental 

( e ) ( p ) 

strains Y, .(At) = Y, , (&t ) + y. (At) and incremental internal state 

(r) ( r )( » ) (r)(p) 

variables (or hidden variables) a 1 ! j (At) = ry t j (Aw + cr t j (At) 
where (a) and (p) represent elastic and plastic components, respec- 
tively. This statement may be given by 


( a ) (p) ( r ) ( a ) (r)(p) 

$(At) = fo[Y u (At), Y u (At), (Y U (At), or tJ (At)] 


1J . . ( a ) ( p ) ( r ) ( • ) ( r ) ( p ) 

a (At) = o[y x j (At), Yij(At), a l i (At), crj j (At) 


( 8 ) 


(9) 


For isothermal conditions, the free energy is the same as the 
internal energy so that 

. . i J. 

f# " Pe - a Y, j 


or for the small time interval At, 


4 (At) = a 1J (At)(Yi](At) + Yj(At)) 


.(e) 


.( p) 


( 10 ) 
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At this point we introduce here the incremental form of free 
energy in a truncated Taylor series expansion, 

1 IJk^ ( e ) ( e ) 1 *1 Jk^ 1 (p) C p) 


P*(At) = 2 E V iJ 


+ 2 E 


l ) 




1 l)d <r)<«> CrHp) (r)(.) ( r )(p) 

+ j ) 5( r ) ( + f *lj K + ) 

r = l 

V ' 11 lJk-2 (r)(a) (r)(p) (•) (p) 

+ P (r) < "k* + »k & )(Y 1 j + Y 1J > 

r=l 

1 J k-fc *1 J k^ 

where E and E represent tensors of elastic and plastic 


(ID 


d 


l 1 k 


9 , 


moduli, respectively; | ( ^ are stiffness constants associated with 

the internal variables. Note that (11) has the form of truncated 
Taylor series expansion only to include quadratic terms. However, 

( e ) ( p ) 

the product term of Y, j and Yj i is missing. This is because the 
coupling of eLastic and plastic strains can he obtained using any 
one of the failure theories and an explicit material kernel relating 

( • ) ( p ) 

the product of Y tJ and Y is nonexisting. 

( r ) 

Lastly, w defined here as the internal variables represent 
time dependent physicochemical properties or simply a viscous be- 
havior which may be expressed as 


= Jexp J V . J 


(r)d - 


( 12 ) 


where r is the time variable and T( p j is the relaxation time. In 
order to facilitate an explicit integration we assume a linear vari- 
ation of Yj j within the time interval At given by 


*.jM- \/*l> 


(13) 


where s is the current time step. Substituting (13) in (12) and 
performing integration we obtain 
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( r 

IT 


s ) 

1 i 


(r)(r), 
= A v <■ 


\ ^ r ^ * / * 
i i ) + n v ( o "i ) 

l j i i 


( r ) . 

+ C Y 


1 i 


( * ) 


(14) 


in which 


( r ) 
A 


exp( 



) 


( r ) 

B 


/r> 

T ( r ) ( D “ 



( r ) ( r) 

C = T (r) (l- D ), 


(r ) 

D 



( r ) v 
(1- A ) 


The derivation of these parameters is given in Appendix 1. 
Rewriting (10) for the current time step as 

. ( S ) • ( e ) . ?$'( 

V. + J- 


\ w?'( 3 ^ s a ) (r). \ 6?’( 8 ) (c)(p) \ -i 

p{— tTT v,/a) + y (,) +- tttt7 v .*, •> + “TTTTTT in <s) l 

* Y i/ s) f> Y i/ s) b ai / 8) fs> 

-J J (.)(V e b> + :. (p) 

v i j 


Yi 


; (8) ) - o 


(15) 


and substituting (14) and (11) into (15) yields 
n 


7*1 


IJk-^ (r)(r)( J . . ( r ) . ( e ) , . ( r )., g . 

) ( A + B Y fcje < •-« + C Y< kX > 


t 1 , .(e) . V* 1 J k £ ( r 1 ' ( r )(p). (»). 

i.))v u (.) +2.S,, < *>h,V 


It) 


r=l 


» , „ . • ( P ) 1 ]|I i.(p),n I 

+ <***> J )Y u (3+ Y u ( = 0 

• (e) 

Since all variations other than Y are not arbitrary we must have 


the relationship 


n 


U, . Uk Jt (,) , V ijki (')(r) , , . (r)*(«), , (r)-(e) . 

<7 3 “ E Y k* * + ^ ?( r ) ^ A Q 'k£ S ‘ 1+ B Y k^ ( ’” 1 5 + C V k 4 <S) ) 

( 16 ) 


r=l 


i 
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*1 p), . 

E Y ( O- a 
i J 


n 

l ]. . • ( p ) T"* 1 J k £ I" ( r ) . 

( ‘ >V U + 2. Sr) L iri 

r 3 ! 


( B ) 


( r ) ( p ) r W/ ( P ) , 

' *tcJ& 




( r ) 


(17) 


Here (16) represents the relationship 


l j b$ 

a = p ~ 77) 

tfij 

which states that the stresses are derivable from the free energy 
functions. It should be noted that, in our specific problem, this 
stress is due to an elastic strain and a law governing the plastic 
strain and stress is needed to obtain the stress due to a total 
strain. The relationship (17) may be considered as the dissipation 
which plays a significant role in heat conduction problems. However, 
for the isothermal conditions as considered in the present study, the 
entire terms of (17) need not be used in the analysis. Only the first 
term will be recovered as we apply a yield criterion in (16). 


II - 4. INELASTIC RESPONSE 

Extensive research has been carried out at Cambridge University 
by Roscoe and his colleagues [19 ] on the subject of the critical state 
soil mechanics. The yield criteria adopted herewere originally pro- 
posed by Roscoe and Burland [19], A plastic tangent matrix in the context 
of the incremental theory of plasticity was derived by the authors 
[22, 23., 2,4 J .a new method of checking conditions of yielding is elaborated 
in Part I of this report. For the purpose of reference we repeat the 
expression for the incremental stress associated with rate-independent 
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elastoplastic behavior, 


1 J 1 Jk^ *1 J k & 

, - <E + E )AV a 


( 18 ) 


A close examination of (16) reveals that a is the total 

stress due to the elastic component of strain and internal variable 
for the current time step. On the other hand, (18) represents an 
incremental stress for a fraction of loading increments with inelas- 
tic strain coupled. It is then immediately clear that if the visco- 

"ki J 

elastic stress as given by (16) is used to calculate E- within 

the time interval and if we proceed with (18) with iterative cycling 
. *1 J k 

for further updating E without participation of the viscous part 
of (16), then at the end of the time interval the total strain reached 
simply reflects the coupling of viscoelastic and plastic properties. 
Thus from (16) and (20), we obtain. 


da U( 


) = 


1 } k-G 
E dy 


k l 


( 3 ) + 


I 

r=l 


1 J k 
H(r> 


^ ( r ) 

( A d 


(r ) 


(s-D 


( r ) . . . ( r ) • , ,« * i J k £ . . 

+ B d^^' 1 > + CdY kfj <«>)+E dY k^ 8> (1-9) 

Note that viscoelastic strain is now associated with the total strain 
as coupling is established. 


II - 5. FINITE ELEMENT EQUATIONS OF MOTION 

The finite element method is widespread in engineering appli- 
cations [11,12]. No elaboration on this method is attempted here. 
In view of (7), (2) and (3) we rewrite (1) as 
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F u dV = 0 

D 


( 20 ) 


n n 

Here the body force pF = F alone is considered merely for sim- 
plicity. The surface traction can easily be added later if needed. 

In the present study we use the plane strain isoparametric 
element with 4 corner nodes. This gives the linear variation of dis- 
placements in the form 


u i = *n u ! < 21 > 

N 

where is the interpolation function and u t is the nodal values 
of displacements u ( (1 = 1,2) and n = 1,2, 3, 4. 

The strain tensor is given by 


'i J 


= 2< u i 


’J 


+ u. 


'!> 


( 22 ) 


Inserting (25) in (26) yields 


Y 


i J 



N 

k 


where 


j( ^ y k 

A ni j = 2 ’(H' n » 1 6 j + Y N >j6i) 


(23) 


(24) 


In view of (21), (23) and (20) we have 


r£w v ^ + Jo u A k N1J dv - J A n <iv}u!! 


(25) 


For all arbitrary values of u k we require the terms inside the bracket 
to vanish, which yields 


••m j n k 
^MN U k A N 1 J 


dv = F 


Nk 


(26) 
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where M MfJ and F Nk are the mass matrix and the force vector, respectively, 


M m n = | PV^V 


i: 

■I: 


F Nk If V .,dv 


(27) 


( 28 ) 


To obtain an incremental form of (26), we take a variation or 


induce a perturbation such that 


M MN dU k 


'f: 


d(T A N 1 j dV d ^Nlt 


(29) 


Introducing the incremental stress (19) into (29) yields 


■M . . £k W ,(«)•£ k (o)£k M , 

M H N dU k (B)+C MAO + ( K «N + Ku > du g (S> 


d F N k (B) 


, <V), V 

+ d V (s) 


(30) 


£ k ( e ) £k ( y ) £ k 

in which C MN , K Mn , and K MN are the viscosity matrix, elastic 
stiffness matrix and plastic stiffness matrix, respectively, 


£k C T* 1 J H n< r ) £ k 

J MN I C )^N»n 

J*7A 


dV 


(*)^K . 


i 

■I 


1 ) a n 


( ^ ) £k | *1 J a n £ k 
*MN “ E A «U A N D n dV 


( V ) 


The pseudo viscous load vector dF Nk is given by 


(V) | ljnn( r )(r), . k 

"/y?<D ^ 

•'v r^l 

ljmn(p) £ k H, 

+ l^ ( r) B A N1J A Hnn dV (du^C.- 


i )' 


(31) 


(32) 


(33) 


(34) 
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The expression (30) is called the finite element equations of 
motion. 


II - 6. SOLUTION PROCEDURE FOR INCREMENTAL 
EQUATIONS OF MOTION 

A solution of (30) can easily be obtained by any scheme of direct 
numerical integration [13]. In this study, a constant acceleration for 
a small time increment is assumed, which gives a recurrence formula for 
displacements , velocities and accelerations in the form, 


f At A t p < e )4k 4* . 

+ 2 C MN + 4 ^ K MN + ' 


dF 


( 8 ) c V ) . 

dF C B ) - O <■ » ) 

ur “'- v Nk 


Nk ‘Nk 


where 


(35) 


. , . H, At s ;,M ( a )^k (y)-^k 

QN ( k s 5 = C HN fdu k (-D + — du k (s-i) 1+ ( k MN + ) 


At a -M, 


{du« ( 8-1 ) + — — du^ 8-1 ^} 


(36) 


M, s •«, v At . At , 

dftjj^ * ^ = du^ s "i ) + — du^^ »-i 5 + — du^ 6 ) 


(37) 


At a M 


du^ 8 ) = du^ 3 " 1 ^ + — — du^ ( s - l ) + —j— <Ju£* 3 5 + Atdu^ °“ 1 ) (38) 


M 

Initially all terms associated with ' s-i) are zero and dii^ 8 ^ in (35) 
can be solved from given initial and boundary conditions. Subsequently, 
du^ 8 ) and du_^ 3 ) are calculated from (38). These responses or his- 
tories are then carried to the next time increment and back to (35). 
However, for the second increment it is necessary to check yield con- 
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ditions and a standard incremental loading method of iteration [14] 
can be applied to each time increment with the total dynamic load bn 
the structure. , 


II - 7. EQUIVALENT DYNAMIC WHEEL LOADS 

Theoretical and experimental studies for the prediction of rigid 
and flexible wheel performance on soil have been reported by various 
authors as mentioned in Introduction. Oriaffeko and Reece [6] 
presented practical procedures in determining radial and tangential 
stresses along the wheel-soil interface. Wong and Reece [7,8] derived 
expressions for sinkage, drawbar pull and torque input based on the 
plate penetration test but with considerations of the important aspects 
of the slip and the actual interaction between the wheels, and spil. 

In the present study the finite element equivalent nodal dynamic 
loadings are determined from the expressions for radial and tangential 
stresses given by Onaf feko and Reece [6] and explicit forms of these stresses 
as elaborated by Wong and Reece may also be used (See Appendix 2 ) . 

In order to compare the dynamic rate -dependent elastoplastic re-.- 
sponses with the results of Perumpral, et al ['12] who neglected the 
effects of inertia and rate-dependency, we consider here the identical 
geometry and material constants. The discretized wheel-soil medium 


is shown in Figure 1. 
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Distributions of radial and tangential stresses on the rim of 
a 49 in. x 6 in. wide wheel on compact sand with 3.1% slip and 41.4% 
slip are shown in Figure 2 and Figure 3, respectively. Equivalent 
nodal loadings as calculated from the tributary area method in Figures 
2 and 3 are shown in Figure 4. It is seen that the area under the curve 
corresponding to the wheel-soil contact area for each finite element 
may be conveniently approximated by the equivalent rectangular block. 

It should be noted that as the slippage increases the vertical down- 
ward loads decrease whereas the horizontal loads increase in the direction 
opposite to the wheel movement. 


II - 8. DEFORMATION AND STRESS FIELDS 


The equations of motion in assembled form for all finite elements 
are solved as described in Section II - 6. In order to compare the 
results for all possible effects, the computer program (Appendix 3 ) was 
written with many optional versions. Various cases studied include 
static analyses for elastic and elas toplastic responses and dynamic 


analyses for elastic, viscoelastic and vis coelas toplastic responses. 


The material constants used are soil modulus E ■ 2000 psi.j Poisson's 

o 

ratio v s = 0.45, angle of internal friction <0 = 36 , density y => 0.05787 
pcij rela xatio n time T = 0.1 sec (r = 1, 2, 3), compression index 

{ *•“) ' ' i 

A. = 0.05, and swelling index h = 0.0001. These constants are chosen 
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to correspond to the compact sand which is used in the equivalent load 
representation as shown in Figures 2, 3, and 4. For dynamic analyses, 
a time increment At = 0.0006 Sec. for viscoelastoplastic response and 
At = 0.0003 Sec. for other responses are used. 

Figure 5 shows these various responses at node No. 31. For 
static analyses, the elastoplastic displacement in the vertical direc- 
tion is slightly larger than the elastic behavior. For dynamic analyses, 
the viscoelastic and viscoelastoplastic responses are considerably 
smaller than elastic and elastoplastic behavior. Once again, effects 
of plasticity result in larger deformations for both viscous and 
nonviscous cases. 

The vector representations of elastoplastic deformations ;for the 
static analysis are shown in Figures 6 and 7. Deformations for 41.4% 
slip are larger than these for 3.1% slip. For the case of dynamic 
analysis (41.4%. slip) the curvilinear transient deformation vectors 
for viscoelastoplastic response are shown in Figure 8. These vectors 
represent the time history from t = 0 to t B 0.6 sec. No doubt that 
the effects of inertia under dynamic loads caused larger deformations 
than under static loads but energy dissipation through the viscous 
behavior retarded, the motion considerably in comparison with the non- 
viscous cases as noted in Figure 5. Deformed shapes for dynamic visco- 
elastoplastic responses at t=0.3 sec. -are shown .in Figure 9-. 



Wheel Contact Length (in.) 

Figure 2:- Radial and Tangential Stress Distribution at the 
Interface for 3.1% Slip on Compact Sand. Ref. [6 
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(T8d) ssaa^s TBT3U38 ubx put? TBjpBH 


Figure 3: Radial and Tangential Stress Distribution at the 

Interface for 41.4% Slip on Compact Sand, Ref [61 


11-70 



NODE NO. 

3.17. SLIP 

41.47. 

SLIP 

f .,(//) 

F„ (« 

v # > 

Fr (#) 

21 

-16.17 

-6.61 

-13.42 

-9.5 

26 

-74.5 

-14.99 

-55.21 

-26.67 

31 

-119.62 

-4.33 

-98.27 

-27.78 

O 4 
JU 

-88.86 

-10. 97 

-89 = 1 

-8,4 

.41 

-27.11 

8.3 

-32.1 

2.85 





i 


Figure 4: Equivalent Nodal Forces as 

determined from Figures 3 and 4. 
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From the deformation fields various stress components are cal- 
culated and the results shown in the form of isobars in Figures 10 
through 21. In the region close to the wheel the major principal 
stresses due to the elastoplastic deformations are smaller than those 
of Ref. [12] as shown in Figures 10 and 11 for 3.1% slip and 41.4% slip, 
respectively. Slightly larger major principal stresses develop at the 
mid-depth for the 3.17. slip. For the case of maximum shear stresses 
(Figures 12 and 13) the present analysis gives larger values than 
Ref. [12] for 3.1% slip, but this trend is reversed for 41.4% slip. In 
general, the maximum shear stresses for the 3.1% slip are larger than 
for the 41.4% slip, the same trend as in the case of major principal 
stresses. Dynamic elastoplastic major principal stresses and maximum 
shear stresses for 3.1% slip at t-0.072 sec., 0.15 sec., H. 228 sec.., .0.3 Tsec, 
0.6 sec. are shown in Figures 14 through 17. Variations of stresses 
with time until maximum stresses are reached are clearly shown. The 
effects of viscosity or rate-dependent plasticity for 3.1% slip at 
t=0.3 sec. and t=0.6 sec. are shown in Figures 18 and 19, respectively. 

The same information for 41.4% slip is given in Figures 20 and 21. It 
is seen that as the slip increases the major principal and maximum shear 
stresses tend to decrease. 

II- 9. CHARACTERIZATION OF SOIL MECHANICS PARAMETERS 

Studies on deformation and stress fields as described in Section 
II-7 indicate that constitutive relationships for the soil behavior sig- 
nificantly influence the response patterns. The mechanics of wheel-soil 
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41.47. Slip) 
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interaction cannot be understood properly if incorrect judgement or 
oversimplification in the theoretical formulation obscures the true 
deformation and stress fields. For this reason, the present study 
was devoted to a new approach in which rate-dependent inelastic be- 
havior coupled with effects of inertia was considered. 

The analysis presented in the previous sections becomes the step- 
stone for characterizing the soil mechanics parameters more realistic- 
ally. Of course, all the results obtained here are based on hypothetical 
material constants. However, if the analytical formulations are cor- 
rect, then the wheel-soil interaction data as observed qualitatively 
and quantitatively may be used to correlate with material constants. 

Such characterization can be achieved by holding some of the material 
parameters constant and comparing the load-deformation data between 
the calculated and observed values. 

Because the present study does not include the dust cloud motion 
behind the lunar rover the observed rooster-tailing cannot be related 
to the material characterization. However, the sinkage of the rover 
wheel together with the vehicle performance data can be used for corre- 
lation with deformation and stress fields as mentioned in the previous 
paragraph. 

II -IQ- CONCLUSIONS 

The main objective of the present study was to introduce a feasible 
constitutive relationship for soil deformation and stress fields under a 
moving wheel. The load transmitted by the moving wheel is dynamic rather 
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than static. The soil is dissipative media in which Inelastic defor- 
mation of the soil is governed by the rate-dependent plasticity or 


viscoeiascOpiaSLiCi-Ly. 


rrf1 J _ 1 J C *.1 - - . r _ __ _ J i J 

ine yie xu oui iduc liicu l y v i ivuowvc auu uui. x auu 


is utilized here for inelastic behavior. The internal variables are 


then introduced to account for rate-dependent viscous behavior. Ef- 
fects of soil inertia are included. Combinations of all of these pro- 
perties result in dynamic analysis of viscoelastoplastic media. 

The numerical results obtained here appear very reasonable. Com- 
parisons with the results of other investigators are made and devia- 
tions are believed to be due to more rigorous treatment of material 
behavior considered in the present study. In order to verify the im- 
pact of the theoretical formulations given here, however, additional 
comparison study through experimental data is needed. 
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APPENDIX 1 

DERIVATION OF INTERNAL (HIDDEN) VARIABLES 


( r ) 


Consider the internal variable <y . . 

1 J 


(r) , X 

«1 j(t) 


= f exp (" ( T ( "^) V U (T)dT 

Jq ' '■ 


(A-l) 


where V 1J ( T ) may be considered to vary linearly within the small 
time interval At, 


Y,,< T > - V,,(t-4t) + - V./t-St)] 


(A-2) 


Substituting (A-2) in (A-l), 
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= exp 
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APPENDIX 2 

CONTACT STRESSES AT WHEEL- SOIL INTERFACE 


The vertical and horizontal forces and torque of a wheel rotating on 
horizontal ground with constant velocity are given by 

9i ^ 9i 


r r 

W = rbf I a(n)co80d0 + I r( 


(0)sin0d0l 


f 01 r 

D = rbf I t(0 )cos0d0 - | CT (0)8in0d0} 

O;, 0 


01 

t(Q ) d0 

in which ct(0) and t(g) are the average radial and tangential stress across 

I 

the wheel width of b(Fig. 2-1). 

The location of the point of the maximum radial stress may be expressed 


-1 


as 

9m “ (^1 + C 3 j ) 0 ! 

where 1 is the slip (%) defined by 


l 


(1- — :)ioo 


and C x and C^, are the constants [14-17] given in Table 1. 


9 
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TABLE 1 


Soil 

Angle of 
Interna L 
Friction 
(deg. ) 

Soil 

Cohe- 

sion 

REBSSj 


C 8 

K i 

K a 

n 

Shear 
Deforma- 
tion Mod- 
ulus (in) 

Compact sand 

33.3 

0.10 

0.0575 

0.413 

0.32 

20 

2.5 

0.47 

mm 

Loose sand 

31.1 

0.12 

0.048 

0.18 

0.32 

0 

2 

1.15 

1.5 

Sand 

36.0 

0.10 

0.0617 

0,285 

0.32 

■ 

n 

B 

- 

Dry sand 

24.0 

- 

- 

0.38 

0.41 

■ 

■ 

- 

- 


In the region between e and q h or the front region, the radial stress 
is given by [18] 


a l (0) = (K 1 +K n b)(-) (cos 0 - cos 9 ^ 

where the constants , and n are shown in Table 1. The radial stress 

acting in the rear region is of the form 

r n / (Cj +C 0 j ) \ n 

ct s ( 0) = (Kj+Kgb) (t) [cosfQ l -0| — — ] -cob 0 1 ] 

\ C i +C 3l / 

The shear stress around the rim is given by [14,15], 


t(0) - (C+a (9) tan 
where C is the cohesion, and 


|c+a(9)tancp^ ^1-e^ 


P 


( 9 1 - 0 )-(l-i )(sin 0 X -sin 0 )} 
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In the above expressions 0 is still not known but can be determined 
from the expression of the vertical force W, 


where , 


>0 i fe* 

W = rb{ I a-t (e)cos 0d0 +1 <j a (8)cos 0dQ 

0 m Jo 


*0i 

+ 1 t x (0>sin 0d0 +| t 3 (e)sin e<*8 


*/ 

-/ 0 ( 


Tl 


/ 

(0) = ^C+ c r 1 (e)tan <jsj|l-ePj 


T 3 (e) - ^c+ ct 3 (0 ) tan |l-e P j 


If the magnitude of W is given then the above integration may be carried 
out by the Simpson's rule and 0 X is solved in terms of known values. 

With the value of 0 X known, we can then calculate the radial and 
tangential stresses. 

Finally, the wheel sinkage Zq is determined from 


z 0 =■ (1- cos 6x>r 
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OF = TPl - NFJ - uFK * (l.-CAPA/RAPD) 

00006300 

ouu 


ASU = SR.b * PP * IPO-PP) 

0000b400 

ouu 


ETA = SQkTITTJ) / HP 

00006500 

UUU 


ViRXTt (b*o2u) VUI0K.T7J* ASGUDF . NEL *SM *ETA 

00006600 

ouu 

62U 

F0R«.1AT(5A* 'VOIbR=* *F10.5* ' 3J=‘*E12.5** ASQS* *E12.5* ' nF=‘* 

00006700 

IJOU 

] 

1 fc, 1 2 • 5 * 17* • R,= *.tl2.5.» ETA=» *E12.5) 

00006800 

ouu 


IF (OR.LT.O.) 00 IU 80U 

00006900 

OUU 


no 3U0 I = 1*3 

00007000 

UUU 


HR l I ) = u. 

00007100 

UUU 


CR < 1 ) = u. 

00007200 

UOU 


UO 3U0 J = 1*3 

00007300 

uou 


PR ( I ) = bRll) + UU,J)*AR(J) 

00007400 

00 u 

.■500 

CR(I) = IRII) + A«(J)*U<U*1) 

onoo75on 

OUU 


DEN - 0 . 

00007600 

uou 


no 310 i = 1*3 

00007700 

uou 

31 U 

CtlM = DEN + AR<I)*bH(IJ' 

00007800 

UOU 


DEN = DEN + bH , ... 

00007900 

UUU 


DO 320 I = 1.3 

0OU08000 

uou 


DO 320 J = 1*3 

00008100 

OCU 

32U 

l’.P (NtL » I * J) = 0(1*J) r Bk(II*CR(J) /DEN. 

00008200 

uou 

0 


00008300 

uou 

HU 0 

CONTINUE 

000084'00 

uou 


no iuo i = i*4 , 

00008500 

uou 


co ion j = 1*4 

00008600 

uou 


ST IFF ( 1 . J) r DP(NtL*l*l)*TYPEA<I * J) +DP (NEL* 3 * 3) *TYPEC ( I **J ) 

00008700 

uou 


STIFF < I , J+4 ) = np(NEL.1.2>*TYPEB(I*J)+DPlNEL*3*3)*TYPEB(D*I) 

00008800 

UUU 


ST IFF ( J+4 . 1 ) = ST1FFU.J+4) 

00008900 

UOU 

100 

STIFF ( 1 + 4 . J+4 ) = UF<NEL*2*2) *TYPEC { I * J) +DP (NEL. 3 *3)*TYPEA(I*J) 

00009000 

uou 

L 

, 

00009100 

uou 


00210 I = 1.8 

00009200 

UUU 


II = KK ( NEL * I ) 

00009300 

uou 


DO 210 J = 1.8 

00009400 

uou 


JJ = KK(NEL.J) 

00009500 

uou 


IFII1.EQ.0 »0R .JJ.LU.O ) GO TO 210 

00009600 

uou 


IFII1.LT.JJ) GO TO 210 

00009700 

uou 


IFdl.GT.IHBI ) GO 10 214 

00009800 

'uou 


L = OJ + UI-1) *11/2 

00009900 

uou 


GO TO 215 

onoioooo 

uoo 

214 

L = JJ + LT + (U-IHB) * IHBl 

00010100 

uou 

215 

XKL(L) = XKL(L) + STIFFU.J) 

00010200 

uoo 

210 

CONTINUE 

00010300 

uou 

900 

CONTINUE 

00010400 


UPU 

0 



00010500 

UPO 



PE IUrK 

00010600 

LI () u 



END 

1)0010700 

000 

la! h L 1 9 

S1H l,Hi A* IPEI..S1FR, ► , 1635601 324 10 


UPO 



SUtO'VI i 1 r.r S 1 F H 

onoooioo 

00 u 



PAKi.'.t [ro i,F l=15u,NFLS=150*MX= 5000 »NF=NFT*2 

OOunOBOO 

Ouo 



CO,.' *,r / i.Lm / 1 ITl-t( 20 ) »lNOr>E»NELEM*NAPC,NF)C»MACT#iSw(b) 

00000300 

uoo 



Uii-k.: /uL.W aKHMX) »P(NF) #IMAX»IhR»IH8I»LT»LA<;l ,NFRFt 

00000400 

000 



C (I|.:'UF /oLi\ 1 '/ 0F» llMtl.S » 3 ) #STRAIB(NELS»3) , ALPHA » BFT A * GAMMA > PtL r » XK * 

1)0000500 

000 



* A 1 , M“0»t . Sr!, SMS, vuIUI »CAPA,RAMD*BET 

QPU00600 

uoo 



COMMOt /uLi\7/ UPdiNF) #U0UB(NF> »FOV(NF> »UDD(NF) »Pd(NF) »L)ISP(NF) , 

00000700 

uoo 



* OUO 1 ( NF ) 

OOOOOBOO 

uoo 



COMMON -/dLKft/ UP UitLS * 3 » 3 ) » SIGB (NELS » 3 > 

00000900 

uoo 



COMMON /UYN/ Cb AK l MX ) * XM(MX ) # AA (MX ) 

onooiooo 

uoo 



NTlME-1 

onooiioo 

000 



NDtLT = 1SW ( 1 ) 

00001200 

000 



N = PACT 

00001300 

uoo 



NPRN f = ISw(j) 

00001400 

000 



IKOUNT = 2S0 

00001500 

oou 



KOUNT = IKOUNT 

00001600 

uoo 



1PLT = 1 

00001700 

uoo 



CALL .ZERU(SIGB,NElS,3> 

00001800 

uoo 

c 



00001900 

uoo 



DO 2 I : l.LAST 

00002000 

000 


2 

A A ( I ) = XM ( I ) + LdAP(I) * UFLT / 2. + XKL<I> * PELT * DELT / 4. 

00002100 

uoo 



DO ISO I = 1 » NFRtc. 

00002200 

uoo 


ISO 

AA ( I+LAST ) = PU) 

00002300 

uoo 

c 



00002400 

uoo 



CALL FAC10FMAA. IHU* IHBl.LT.LASTfNFHEt) 

00002500 

uoo 



CALL SOLTN (AA, IHd» IHB1 * LT » LAST » NFFiEE ) 

00002600 

unu 

0 



00002700 

uoo 



DO 3 1=1, N 

00002800 

uuu 



UDb(l>=0. 

00002900 

uoo 



UDUB ( I ) =U . 

00003000 

uoo 



FOv ( 1 ) =0 • 

00003100 

uoo 



UDU(i) = AA(LAST+1) 

00003200 

uoo 


3 

DBII)=0. 

0O003300 

uoo 



DO 4 1=1, N 

00003400 

uoo 



DlSP ( I ) = Dd ( 1 ) + UELT *UOb ( I ) +(UDDd(I> +UDD < I ) ) *DELT**2/4 . 

00003500 

upo 


4 

DDOT ( I ) = UUP(I) + (UUDd.(l) +UDPU) )*DELT /2. 

00003600 

ouo 



WHITt (6, 100) NTIMt 

00003700 

UPO 


loo 

FORMAT (25X, 'DISPLACEMENTS FOR TIME INCREMENT » r 15 • 5X ,* LINEAR » > 

00003800 

uoo 



fcPiTt (6,6947) (OlSP(I) ,1=1,N) 

00003900 

uoo 


6997 

FORMAT (8E1S.6) 

00004000 

UPO 


101 

FORMAT (BtlS.o) 

00004100 

Ouo 



toMTE(6»102> 

00004200 

uoo 


102 

FORMAT ( ’VELOCITY' ) 

00004300 

uoo 



VvPiTt (6,101) (UPOl (I) ,1 = 1,N) 

00004400 

uuo 

o 



00004500 

oou 

0 


INITIALISE Ot> , STHAlIT 

00004600 

IJOO 

L 



00004700 

ouo 



CALL ?EK0(STHAlB,NtLS*3) 

00004800 

uoo 



CALL 7LK0 ( o*-> »NLL5 , J ) 

00004900 

000 

0 



00005000 

uoo 



CALL ST An ( r,TlMt > 

00005100 

Uuo 

c 



00005200 

uOu 

L 


NOw READY FOR l\Ew INCREMENTS.. 

00005300 


uoo 

0 



00006400 

uou 


10 

CONTINUE 

00005500 

000 

c 



00005600 

000 



Jd = IHR1 

00005700 

000 



DO 200 1 = 1 fNHPtu. 

onoo5flno 

ooo 



JJ = JJ + 1 

00005900 

000 



SUN = 0. 

00006000 

000 



SUN = 0. 

00006100 

ooo 



IF(l.GE.lHb) GO TU 180 

00006200 

000 



L = (I+l)*i / 2 - i 

00006300 

000 



DO 170 J = 1 » do 

00006400 

ooo 



IF(sJ.LE.l) L = L + 1 

00006500 

000 



IF<d.GT.i.ANU.d.LE.IHB) L = L + d - 1 

00006600 

ooo 



IF ( J.GT .lHb) L = L + IHB1 

00006700 

ooo 



SUN = SUN + XMIL)*0B(J) 

OOOO 6 BOO 

ooo 


170 

SUN = SUM + LRaR IL) * (UDB < J) +UDDR t J) *DELT/2 . ) ♦ XKL (L> * ( DB ( J ) +DELT*00006900 

uoo 


lUDb ( u ) +UUDB ( J ) *DElT*0ELT/4 . ) 

00007000 

ooo 



GO TO 199 

00007100 

OOO 


180 

II = I - IHB + 1 

00007200 

ooo 



L = LT + (11-1) * IHB 

00007300 

ooo 



IF( Jo.GT .NFREE) dd = NFRtE 

00007400 

uoo 



DO 1.90 d = II. dJ 

00007500 

uoo 



IF(J.LE.l) L = L + 1 

00007600 

ooo 



IF(J.GT .1) L = L + IHBI 

00007700 

ooo 



SUN = SUN + XM(L)*UB(J) 

00007600 

ooo 


190 

SUM = SUN + CBAR(l)*( UUP ( d ) +UDDR ( J ) *DELT/2 • ) + XKL (L ) * ( Db (d ) +DELT *00007900 

uoo 


lUDb (o ) +UUDB ( J) *DtLT *UELT/4» ) 

00008000 

ooo 


199 

XM (LAST + 1 ) = SUN 

00008100 

UOO 



CHAKU AST + 1 ) = SUN 

00008200 

000 


200 

CONTINUE 

00008300 

uoo 

c 



00008400 

uuo 



IF (NTIME.GT .NDtLT i GO 10 113 

00008500 

ooo 



DO 12 I = l.NFREt 

00008600 

U(,0 


12 

AA ( I+LAST ) = Fill - FOVII) -XNtLAST+I) - CPAR(LAST+1> 

00008700 

uoo 

c 



ononoflon 

ooo 


300 

CALL SOLTN ( AA . iHb . IHB 1 . LT » LAST . NFKF.L > 

00008900 

UOO 

c 



00009000 

ooo 



DO 13 1=1. N 

00009100 

ooo 



FOV ( I ) =0 • 

00009200 

uoo 



UDD(l) = AA(lAST+1J 

00009300 

uoo 



OISP(I)= Dti ( 1 ) + L/t.1 T *UDb ( 1 ) +UIUrib(l) +UP()(l))*DELT**2/9. 

00009400 

uoo 


13 

ODOT ( I ) = UUB(I) + UlUDbd) +unn<l ) )*DELT /2. 

00009500 

ooo 

c 



00009600 

ooo 



CALL STAN(NTlNt) 

00009700 

ooo 

c 



00009800 

ooo 



IF (NTIME.Nt.NPHNT 1 GO TO 310 

00009900 

ooo 



NPKNT = NPKNT + lbW(3) 

onoloono 

ooo 



WRITE (6. 10b) NTIME 

oooinioo 

onu 


105 

FORMAT ( 1 OX . ‘UISPLaCENENT AT NTIME =’ . I5/5X. ’NODE NO Z-DISPL 

00010200 

ooo 



1 K-DISPL’I 

oooio^on 

ooo 



DO 106 I = 1 . INODE 

00010400 

000 



dd = I * 2 

00010500 

ooo 



II = dd - 1 

00010600 

ooo 


106 

WR 1 TE ( 6 . 620 ) I » OISP < I I ) » UISP ! JJ > 

OOOJ 0700 

000 


620 

FORMAT! 19. 2E 15.6) 

00010800 

ooo 



11 = LAST + 1 

00010900 

000 



dd = LAST + NFREE 

00011000 


1 Ml u 


V*H 1 Tt ( 6 . taOU ) 

onoii ion 

UllU 


WPirt(b.iOl) ( o A 1 1 ) . 1=1 I < JJ ) 

001)11200 

uuu 


H. K i 1 1 < 6 » 0 1 U ) 

001)11300 

UOU 


WRITE (6. 101) (CBMKlI) . 1=1 I » Jo) 

OOU114Q0 

UllU 

60U 

FORMAT (/bX» * fd .FORCE'/) 

onoi lsoo 

uou 

61U 

format ( /i>X . * FP • / ) 

OOU1160O 

Ouu 

31U 

CONTINUE 

00011700 

uou 


GO TO 10 

onoiiaon 

UllU 

c 


onulisoo 

UOU 

113 

WKI Tt ( 6 » 19t>8 ) 

00012000 

uou 

1968 

FOkMAT (//2X» 'END uF ANALTSIS'///) 

OOU12100 

uuu 


RETURN 

onu 12200 

uou 


ENu 

00012300 

Uuu 

wELT » 

SIH NASA*TPF*.STAN» > » 163570132410 


UllU 


SUbRONTlNE STAN(NilMt) 

ooooolon 

ono 


PARAMETER NFr=15U»NELS=lbO.MX= 5000 . NF=NFT*2 

00000200 

UOU 


COMMON /ULKO/ Tnt-£(20)#lNCint»NFLEM»NAPC*NHC»VACT»ISw<5) 

00000300 

uuu 


COMMON /BLK2/ inii,F .2) . IdKL(NFLS*4> .H(NFT> . Z (NFT ) »KK(NELS,a> 

0 0 U 0 0 4 0 0 

uuu 


COMMON /bL*3/ XKLIMX) * F* ( NF ) .IMAX.IFiR. TFIRI *LT .LAST » NFrFE 

00000500 

uou 


COMMON /BLK4/ U(3»3> * ATO (3*3) »A0(3»3) .BD(3.3> * ST IFF (8*6) rCM(«.b) » 

00000600 

uuu 


* Vt(a.b) . Ai (MELS) 

00000700 

UUU 


COMMON /uLNS/ OH (nELS . 3 ) .STRAlR (NFLS.3) » ALPHA . HFT A ► GAMMA » nEL I »XK» 

00000400 

000 


* AT » XNU . E * bM » SMS * I/O I Li I . C ApA ,H AMD . BlT 

oonniiPon 

uou 


COMMON /0LK6/ bT (l’«tl S » o » 3 ) » AKM < NELS * 4 ) .AZMINELS.4) .AOJ(NbLS) 

ononioon 

uou 


COMMON ZuLn 7/ 00b INF ) »UDUH(NF ) >FOV INF ) » HOD (NF ) .OB <NF ) .UlSp ( NF > . 

oouni inn 

UOU 


* liLKJT (NF ) 

00001200 

uou 


COMMOf /uLkF,/ uP<NtLS.3«3) »SI6R<NELS.3) 

00001300 

uou 


COMMON /uYIm/ CBAHIMX) »XN(MX) *AA(MX) 

00001400 

uou 


DIMENSION V(4) *11(4) .VI) (4) ,Un(4> .EPS<3> »FFSD(3> »S1GE(3> *SlGV(3> . 

00001600 

uuu 


* S1GF (3) . bloT ( 3 ) * F V ( 8 ) 

00001600 

uuu 

c 


00001700 

uou 


CALL ZtRO(XM.MX.l) 

ononioon 

uuu 


REm INC 2 

00001900 

uou 


NN = NT I ME / ISfc <o) 

00002000 

uou 


MM = NN * ISV»(J) 

00002100 

uou 


IR1TE = U 

00002200 

uou 


IF (MM. EG. NT IME ) iixlTt = 1 

00002300 

uou 


IF IMP. EO. NI If, El fchiTt (b.oOU) NT IME 

00002400 

uou 


DO 900 NtL = 

00002500 

uou 


00 1U0 I = 1.4 

00002600 

uuu 


IN = IJKl(NEL.1> 

00002700 

uou 


II = (1N-1)*2 ♦ 1 

onoo 2 «on 

uou 


vui = uisp(ii) 

00002900 

IHIU 


U(l) = D1SP ( I I + 1 ) 

00003000 

UOU 


VD(I) = UnOT(Il> 

00003100 

uou 


UU(I) = LOOT ( I i + 1 ) 

00003200 

uou 

10U 

CONTINUE 

00003300 

UOU 


CALL ZERO ( EPS *3.1) 

(10003400 

uou 


CALL ZERO ( EPSD » 3 . A ) 

00003500 

uou 


DO lin I = 1.4 

00003600 

uuu 


EPb(l) = EPSU) + ARM.(nEL»1)*V(I> / AOJ(NEL) 

00003700 

uou 


E PS ( 2 > = EPS 1 2 ) + AZM (NEL » I ) *U ( 1 ) / AOJ(NEL> 

00003400 

Uou 


EPS ( 3 ) = EPS(3)+(aKM(NEL.I)*U<I)+AZM(NEL»I)*V(I)) /AOJ(NEL) 

00003900 

uou 


EPSU(l) = EPSD ( 1 ) + ARM(NEL,I)*VD(I> / AOJ(NEL) 

00004000 

uou 


EPSD (2) = EPbO ( ? ) + AZM (NEL . I ) *UD ( I ) / AOJ(NEL) 

00004100 

UOU 

11U 

FPbll(3) = EPSD(3)+(ARM(N£L»I)*U0(I>+AZM(NEL»1>*VD(I> ) /AOJ (NEL ) 000(14200 

uou 


00 120 1 = 1.3 

00004300 



000 ' 



SlbE(I) = u. 

00004400 

000 



SIfaVtl) = u. 

00004500 

000 



SIbFd) = 0. 

00004800 

000 



DO 119 J = 1*3 

00004700 

000 



SlbE(I) = SIGEd) + 0<l.j)*FPS<J) 

00004800 

000 



SlfaVlI) = blfaV(I) + ATud*J)*EPSO<J> 

00004900 

000 


119 

SlfaF(I) = SIGFd) + bOd*J)*STRAIB(NEL*J) + AD ( I . J ) *bR t NEL * J ) 

()O()05O00 

000 


120 

SlbTlI) = SIGEII) + blbV(I) + SI6F1I) 

00005100 

000 

c 

' 


uoo 05200 

000 

c 


COMPUTE PRINCIPAL STRESStS 

0 n 005300 

000 

c 


COMPRESSION IS PUS11IVE HERE 

00005400 

000 

c 



00005500 

000 



SOK = 1 1 (S1GT (l)-SAGT (2) )/2. )**2 + SIGT<3)**2) **.5 

00005800 

000 



SGI =-(SlGTU> + SlGT (2 ) ) / 2. 

00005700 

000 



EPS IT) = So I + SuR 

08005800 

000 



EPS (2) = SGI - SUR 

00005900 

000 



EPS (3) = SUR 

OOOObOOO 

000 



IF (MM .EQ.N1 IME ) VvblTE 1 6 * bl 0 ) NEL * S1GE * SIGV , SIGT * EPS 

ijouom on 

000 

c 



onuobPnn 

000 



IF1ISW12) .EO.l) ChLL PbTlFMNEL*SIfaF*IPLST*IRTTF) 

(10008300 

oou 



IF(ISW(2).tO.O) GU TO 90U 

onuof'4nn 

000 

c 



onor,b5oo 

000 

c 


CIJMPUIE Vise. t-ORCE FOR EACH ELtMtNT 

oooObGon 

000 

c 



onoiih7rio 

000 



DO 130 1 = 1*8 

00005800 

000 



F V < I > =0. 

i)8ooh Q on 

000 



DO 130 J : 1*3 

00007000 

000 


130 

F V ( I ) = F V ( I ) + bllNEL*I*J) ♦ SIGF(J) 

000071 00 

000 

c 


NUh ASStMbLt in GLObAL FORM 

00007200 

000 



DO 140 1=1*8 

00007300 

000 



II = KK(NEL*i) 

00007400 

000 



IF ( 1 1 . bti • 0 ) bO TU 14 U 

00007500 

000 



FOV(II) = FOVlllJ + FV(I) 

00007800 

000 



IF (ISW12) .nE.I.Or.IPLSI .NE. l) GO TO 140 

()8(J077()0 

ooo 



DO 1 39 J = 1 *8 

08007800 

000 



JJ = KK ( NEL * J ) 

1)00079(10 

000 



IF ( Jo .to . 0 ) bO TO 139 

0 8 0 0 ri 0 0 0 

000 



IF1I1.LT .JJ ) Gb lu 139 

ooooHion 

uoo 



IFdl.fal.lHPl) GO 10 137 

08(108200 

ono 



L = oJ+(lI-l)*lI/<: 

00(108300 

oou 



GO 10 13b 

00008400 

000 


137 

L = oJ + LI + 1 1 1-lHb ) *IHH1 

0OUO850O 

ooo 


138 

XML) = aM(L) + SliFF(l*j) 

00008600 

oou 


139 

COhl II ! OE 

00008700 

uoo 


14U 

CONTINUE 

00008800 

ooo 

c 



00008900 

0(10 

l 


UPDATE (GO) Anu (STKA1B). (ATRA1N RATE) 

00009000 

ooo 

L 



00(109100 

ooo 



DO 150 I = 1*3 

00009200 

ooo 



Slbrt U’lEL * I ) = -blut(l) 

00009300 

uoo 



GHINtL.l) = ALPHA *OH (NEL *1 )+BETA*S1RAIB< NEL * I ) +GAMMA*EPSD ( I ) 

(10009400 

0,00 


lbO 

SlRAiBlNcL* I ) = tPSDd) 

00009500 

oou 

L 



00009800 

uoo 


900 

CONT 1NUE 

00009700 

oou 

c 



00009.800 

uoo 

c 


OPDATt OISPL . * vtL . * AND ACCEl.L. 

00009900 

uoo 

c 



oooioooo 


UOU 
U(.U 
U Cl u 

OCIU 

U (. u 
(J (: U 
Ul'lO 

uou 

uou 

Ul'l) 

U()U 

uou 

UuO 

uou 

uou 

uou 

uou 

uou 

uou 

uou 

uou 

uno 

uou 

uuu 

uou 

uou 

uou 

uou 

uou 

uou 

uou 

uou 

uou 

oou 

oou 

uou 

oou 

oou 

uou 

ouu 

oou 

000 

000 

000 

ooo 

000 

000 

ooo 

uou 

uou 

^lou 

uou 

000 

ooo 

ooo 

oou 

uoo 


lhU I'-liVt = NTU't + 1 
1(0 sir i = 1 1 mac t 
PHI) = UlSP(I) 

UDbU) = nuOT(l) 

91 u UUulilJ) = Uf:u(l) 

6UU FOnRM (///10X» ’NfihK = • * 1 5/2X » * ELEN' * * T1 0 , * ELAST IC STrESS'»T40» 

1' VILLOUS SIKtSb* * l 7 Of *101 .STRESS' » T 3 U 0 » •PRINCIPAL STRESS'/) 

61 U FORMAT ( lb» 12tlU.3) 

L 

620 FORMAT (//10X » 'PRINCIPAL STRESS. MURE =«»I5/) 

630 FORMAT ( 1 10 » 3E13 • 5 ) 
l 

RETURN 

EMO 

laLLI »SIH NASA*TPF$.PSTIFP, t >163602132410 

SUBROUTINE PSTIFK INF L »SIGT* IPLST # IKITE > 

PARAMETER NFT=16u»NELS=lbO»MX= 500U.NF=NFT*2 

COMMON /t)LKO/ UTLt (20 ) . INODE .NELEM r N APC . NBC * MATT f I Sw < b ) 

COMMON /h)l-K 1 / nv<4)*H(4)»AR(4)*BR(4)»CR(4),A?(4),HZ(4)rCZ(4>» 

* f'N ( 4 ) .Ch<4) fDn( 4J » T YPt A ( 4 > 4 ) »TYPEU(4.4) .TYPtC(u,4) » TYPED (4 » 4) . 

* A(j > BO • Co i 1 < ( J4 » Rv- r l.C * MEL 

COMMON / Ll. i / b(j. J) »A1D<3*3> ,AD(3»3) iUDn,3) »ST1FF<8»8) fCM(8»8) r 

* VL ( o • fa ) » ' - ' lELS ) 

COMMON /Bl.rO / ult(„Ll.S»j) *S1RAIR(NELS»3) > ALPHA . RFT A » GAMMA • DELT . XK t 

* A r r ANU»L»Sr *Sl«S, \,U1DI »C«PA »RAMD#BtT 
COMMON /riLrvr/ UP(iMtl.SOO) » SIGB < NELS » 3 > 

DIMENSION SIGTC3) 

C 

c compression is positive 

c 

CALL ZERo<STlFR.a»a) 

READ ( 2 ) f YpEA * TYPcb » TYPEC t TYPED 
1PLST = -1 

IF (NEL.Lt . ISW ( b ) ) 00 TO UOU 
XNUS = XNU * XNU 
V A - XNU + 1. 

VB = 1 .+XNUS-XNU 

VC = 2 . * ( XNUS-XNU ) - 1. 

DSIGZ = -S1GT ( 1 ) -SiGB (NEC 1 1 ) 

DS1GR = -SIGT(2)-S1GB <NEL 1 2 ) 

DTAUZ = -SlGT<3)-SiGB<NEL»3) 

SIGZ = -SIGT(l) + US1GZ/2. 

SIGR = -SIGT (2) + USIGR/2. 

TAUZ = -SIGT (3) + UTAUZ/2. 

DELP = VA*(DSIGZ+uSIGR) / 3. 

PP = VA * (SlGZ+SiGR) / 3. 

PS6l = PP * PP 

SZZ = 2 • *VB*SIGZ + VC*SIOR 
SRR = 2 • *VB*SIGR + VC*S1GZ 
C -/R = 6.*TAUZ 

io = (VB*SIGZ*SIGZ+VB*SIGR*SlGR+VC*SIGZ*SlGR>/3. + TaUZ*TAUZ 
Vi J = 3 . *TJ 
i 1 1- = l.+TTJ/<PSQ*SMS> 

POn =: l.-CAPA/RAMU 
PO = PP*ETM**POW 
CALL tREAA(NEL* AREA) 

RA| 10 = AREA / A 1 INEL ) 
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Vll< = VOiDh + 1. 
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AA = SMS * ( <; ••♦PR-PO ) / 3. 

00004500 
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00004700 
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(inoopi on 

UOO 



IU' (3) = SZR 
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00006400 
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DO 1 = 1*3 

1)0006600 
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111- i = UFo + MR l I ) * STrA iH (NLL * 1 ) 
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0(H) 
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oil (I 
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IF (NP'.LT.O. ) 60 10 UOU 
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DO 300 I = 1*3 
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PPII) = u. 
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CR ( 1 ) = u. 
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no 3uP J = 1*3 
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E'Fill) = uHlI) + UU.J)*Ar<u) 

00007300 
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C R ( I ) = Him + Ai\ ( J ) *U ( 0 * 1 ) 

00007400 
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DO 310 1 = J .3 

00007600 
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DEN = DEN + ARlT)*bR(I) 
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DEN = DEn + bF> 
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DO 320 I = 1*3 
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DO 320 J = 1.3 
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DP ( MEL . I * J ) = -RRli) * Ck(J) / DEN 
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00008500 

uoo 



no ion j = i*4 

00008600 

uoo 
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ST IFF ( 1+4. J+4 ) = L/P (NEL *2 * £ ) *TYPEC (I * J ) +L)P ( NEL * 3 * 3 ) *TYPEA ( I . J) 
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END 
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DIMENSION XKll) 
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N = PiFREt 

00000400 
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IHol = IHB I 
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00 * 1=1 iM 
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IF (I. OT.iH.il ) on iu 2 
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K = 1 

oooouson 
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M=K+(I-l)*i/2 

onoooooo 
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K=l-iUHl 
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F=K+LT+( i-lHH)*Ihttl 
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b=i+lUbl 
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unooi non 
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JJ=I+IhHl 

ooooi son 
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Go To b 
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00001700 
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H = U • U 
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l A=1 -1 

uoooi^oo 
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LI =1 + 1 

00002000 

Uoo 



IF(LA.Eb.l)) oO To o 

00002100 
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no / L=K.La 

1)0002200 
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IF (L.01 . 1 HtJ 1 ) GO lo SO 

00002300 

00 0 



J = (L + i ) *L/2 

00002400 

000 



GO TO bl 

1)0002600 

000 


bO 

J = LT + IHH* ( L-ii Ifa] ) 

00002600 

Ouu 


bl 

A = XK (M) 

00002700 

uoo 



H - H+A*A*XM«) 

00002600 

uoo 


7 

Mrh + 1 

001)02900 

00 0 


fa 

A = XK. 1 M ) 

00003000 
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XK (M)=A-t) 

00003100 

ooo 



IF(I.FQ.N) GO 10 o 

00003200 

000 



no 9 J=Lb . JO 

uno03300 

0(10 



SUM=O.U 

00 003400 

uuu 



IF(J.GT.iHbl ) GO 10 10 

00003600 

ooo 



K = 1 

00003600 

uoo 



MM=K + ( J-l ) * J/2 

00003700 

UUU 



GO 10 11 

00003BOO 

ooo 


10 

K=u— I HB 1 

00003900 

uou 



MM=K+l T+ ( 0- IHP ) * Xnb 1 

00004QOO 

uoo 


11 

IF(LA.EQ.O) GO TO 9 

00004100 

uoo 



IF (K.G1 .LA) GO TO 9 

00004200 

ooo 



00 12 JA=K.LA 

00004300 

uou 



L=m-1+JA 

00004400 

uou 



IF ( JA.GT. IHD1) GO 10 13 

00004500 

ooo 



Ll= ( OA+l ) * JA/2 

00004600 

ooo 



GO 10 14 

00004700 
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13 

Ll=LT+lHb* l JA-IHrci ) 

00004fl00 
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14 

SUM=SUM+XK( MM )*XK(L)*Xt\(Ll) 

00004900 
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12 

MM=MM + 1 

00005000 

000 


s 

XK (MM)=(XKlMM)-SUn)/XK(M) 

00005100 
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a 

CONTINUE 
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RETURN 

00005300 

000 



END 
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THIS PORTION OF SubPOUT INE PERFORMS FORWARD-SUBSTITUTION 
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DIMENSION XK ( 1 ) 
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N = NFREt 

00000500 
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iHbl = IHB I 

00000600 
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NF = LAST + 1 
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DO 1 K = 2»M 
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000 
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ooooioon 
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IF (K .61 . lHbl ) 00 10 2 

00001100 

000 



M=U 

00001200 

000 



MM=K-1 

00001300 

000 



Ml=MM*K/2 

00001400 

000 



60 TO 3 

00001500 
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M=K-1HB 

00001600 

000 



N,M=IHR1 

00001700 

000 



M1=M* IHB1+LT 

uoooikoo 
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SUMrO.O 

00001900 
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DO 4 L=1»MM 

00002000 

000 



LL=L+M 

00002100 

000 



ju=ll+mi 

00002200 

000 



LL=LL+NF-1 

00002300 

000 


4 

SUM=SIJM+XK t Jo) *XK ILL ) 

00002400 
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1 

XK(LL+1)=XK(LL+1)-SUM 

00002500 

ooo 



J - nF+N-1 

00002600 

000 

L 
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00002700 
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MF=NpTN-i 

00002800 

ooo 



XK (NF ) = XI\(NF)/XK (LAST ) 

0O002900 

OUO 



DO 5 K=2#N 

00003000 

ooo 



l=n-k+i 

00003100 

ooo 



IFIL.GT.IHBJ ) 60 106 

00003200 

ooo 



I=L+(L-l)*L/2 

00003300 

000 



60 TO 7 

00003400 

000 


6 

I=L+(L-IHF3)*IHbl+LT 

00003500 

ooo 


7 

IRrN-IHB 

00003600 
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IF(L.GT.IR) GO TO b 

00003700 

oou 



J=lHbl 

00003800 

ooo 



GO TO 9 

00003900 

ooo 


6 

J“h-1 

00004000 

ooo 


9 

SUM=U .0 

00004100 

ooo 



DO lu M=1»J 

00004200 

ooo 



NM=L+U 

00004300 

ooo 



IF(MM.GT.IHHl) GO TO 11 

00004400 

ooo 



NN=L+(MM-l)*MM/2 

00004500 

ooo 



GO TO 12 

00004600 

ooo 


11 

NN=L+ (MM-IHB) *IHbi+LT 

00004700 

oou 


12 

MM = NF -N+MM 

00004800 

ooo 


10 

SUMrSUM-MK (NN)*XK (MM) 

00004900 

uou 



mm=nf-n+l 

00005000 

oou 


5 

XK ( MM ) -Xt\ ( MM ) /XK ( 1 ) -SUM 

00005100 

oou 



RETURN 

00005200 

ooo 



ENU 
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SUBROUTINE SETUP 

onoooloo 
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COMMON /oLM/ W<4) rH(4> »AR(4) »RR(4) »CR(4 ) ,aZ( 4> ,BZ(4) tCZ<4> » 

00000200 

ooo 



* BN ( 4 ) r CN ( 4 ) > ON ( 4 ) rTYP£A(4»4) »TYPEB(4.4) »TYPEC(4»4) *TYPEU(4»4 ) > 

00000300 

uou 



* AOrbOrCO* iCr JCfKCfLCrNEL 

00000400 
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W(l) = 0.3476548 
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ooo 



W( 2) = 0.6521462 

00000600 

oou 



W ( 3 ) = In ( 2 ) 

00000700 

ooo 



Vw ( 4 ) = Will 

OOOOOROO 
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H ( 1 ) r 0.8611363 
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M2) = 0.3399810 
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H ( 3 1 r -nU> 
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uoo 


H(4) r -H ( 1 ) 
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HM1) = -1. 

unuo l 3on 
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PN ( 2 1 = 1. 

00001400 

UOU 


HN13I = i. 

onumsno 

uou 


B N ( 4 1 = -1. 

uoooiaoo 

uuu 


CNU) = -1. 

00001700 

UOU 


CN(2) = -1. 

ononi«on 

uou 


CN ( 3 1 = 1. 

ononionn 

OOU 


CM (4) - 1. 

00002000 

uoo 


LiM(l) = i. 

00002100 

uog 


CN<2) = -1. 

0011112200 

uoo 


Dig (3) = i. 

00002300 

uou 


DM 14) = -1. 

00002400 
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TW0P1 = 1. 

onunpson 

uou 


RE 1 URN 

onoo2Bnn 

uou 
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0000271)0 

uuu 

NEL 1 t 
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COMMON /dLx, 1 / w<4 1 »H(4) #AR(4) »PR(4) »CR(U) ,A7(4) ,B7(4) »C7(4) » 

00000200 

uuu 


* H M ( 4 ) »CW ( 4 ) »Dn I 4 1 « T YPt A l 4 » 4 1 »TYPEB(4»4) »TYPtC(4»4) *TYPEU(4#4) » 

00000300 
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* AO f BO # CU r 1C # JC # KO # LC # MF L 

00000400 

UOU 


IPT = 4 

ononosoo 

OOU 


AM = U. 

onunoBon 

UOU 


00 lUO I = 1 » 1HT 

onoon7on 

uoo 


X = Ml) 

ononoRnn 

uou 


DO lun U = 1 » IHT 

on onoRon 

uou 


Y = H(J) 

nnumonn 

uou 


AA = AA + W(l> * *(J) * HKiXiYI 

uooni l on 

uou 
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C0MT1MUE 

00001200 

uou 


PE IURN 

00001300 

uuu 
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00001400 
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FUNCTION FlK.X.Y) 

ononol on 

uou 


COMMON /bLM / w (4) #H(4) #AR(4) #HR<4> #CR(4) ,AZ(4> ,BZ(4) #CZU) # 

00000200 

uuu 


* I'N (4 ) # Cl\l ( 4 ) » DIM (41 »TYPtA(4f 4) #TYPEB<4#4> #TYPEC(U»4) »TYF’FU<4»4) * 

1)0000300 
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* AO t bO » Co * 1 C * JO » KC # LC # NEL 

1)0000400 

uou 


COMMON /bLKR/ M * N 

ononosoo 
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FC = AO + bO*X + 00* Y 

onoooBoo 

uou 
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00000700 
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GO TO (1.0#20»30#4U*b0»6u) * K 

oooooBon 

uou 

10 

CONTINUE 

ooonoRon 
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AM AIM = AR (m) *AK <|M) 

ooooinon 

uou 


RMAN = Hr(M)*AK<N) *■ UK<N)*AR(M) 

onooi ion 

OOU 


AMCN r AR ( M ) +CK (NX + AK<N)*CH<N) 

00001200 

uou 


PMBN = PR ( M ) *BK ( N ) 

00001300 

uou 


BMCN = HR ( M ) *CR ( im ) + CR(M)*HR(N> 

00001400 
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CMCN = Ck<M)*CR(W) 

onooisoo 
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GO TO 1UU 

onooifsnn 

uuu 
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continue: 

00001700 
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AMAIM = Ar(M)+AZ<M 

OoooiROn 

uou 


BMAn = UZ(N)*Ar(M,) + HR<M)*AZ<N) 

ononiRon 

uuu 


AMOM r Ar(m)*CMM + CR(M)*AZ<N> 

00002000 

■ uou 


EJMHN = BR ( M ) *H2 ( N ) 

00002100 

uou 


PMON = BR(m) * CZtN) + CR(M) * BZ<N) 

00002200 

uuu 


CMUM = CzOmUCKIm) 

00002300 

uuu 


GO TO lOu 

00002400 

Uou 
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CONTINUE 
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AKmN = AZ<M)*AZ(N) 

BR.mN r bZ<M)*AZ(N) + B^<N)*AZ<R> 

ARCN = AZ<M) *Cz<h) 4 ( tM ) *Ci: < N"- ) 

PNbN ^ BZ (N ) *BZ (N ) 

BRCN r BZ(M)*CZ(N) 4 CZ <M > *HZ ( N > 

O'UJ = CZ(M*Cz<N) 

GO TO illO 
40 CONTINUE 

XY = X * Y 

F A - ] • 4 bN<M)*X 4 CN < N 1 ) * Y 4 RN(M)*X Y 
OH = 1. + bN (N ) *X + CN IN ) *Y 4 DN<N)*XY 
F = FA*Fb*FC/lZ0. 

RE IUKK 

SO F = IAB(M) 4 BK(M)*X 4 Ch ( M ) * Y ) / FC 
RFHlRh 

00 F = (A2<M) + BZ(N)*X 4 Cz(fo)*Y) /FO 
PfltlHF 

100 FA = A MAN + rtMAN* A 4 AMCN*Y 4 BMbN*X*X 4 Bn*CN*X*Y 4 CMCN*Y*Y 
F = F A/F O* . 1 25 
KF 1IJKN 
F MU 
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SUdROI IT 1 NE t/lAOML UK f AHF » NKRtE » I H« » I HR I »LT,LA4T) 

PIMENSION XK (1) .ABF ( 1 ) 

0 

no 200 1 = lfLAST 
2UU XK ( I ) = U. 

0 

no 100 J = 1 * Nl-Rtc. 

IFlJ.GT.lHbl) GO lo 10b 

I = lv)+l) * 0 / 2 

GO TO 109 

10b L - IT 4 IHB * (o - IHbl ) 

109 XKIL) = i. 

1UU CONTlfUE 
O 

no 110 I = ) iNl-Ril 
11U XK(LwST41) = AFF(i) 

0 

RETURN 

ENU 
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COlsMOf' /bLk 2/ lOlbF.2) f IsJKL(NELS>4) rR(NFT) r 7 ( NFT ) »KK (MELS , b ) 
COhMOT /dLK7/ ODbUJF ) t ONUB 1 NF > »FOV(NF > tUUD (NF ) »Db (NF ) »UISP (MF ) t 
* nU01(NF) 

DIMENSION r<R(4) »Zz_<4) 

no iuo l = ) »4 

IN = IJKL(NFL.i) 

JU = IN * 2 

II = JJ - 1 

PR ( I ) = K<lN) 4 uibP(JJ) 

100 ZZ (I) = Z<1N) + UibP(Il) 

A I = (RR(2)-RRll ) )*(Z7(4)-ZZ<) > >-<RP(4)-RR< 1 ) )*<ZZ(2)— zz(l ) ) 

A J = (RR(3)-KRl?) >*IZ7AU)-ZH?) )-(HP(U)-RR{?) >*(ZZ(3)-ZZ(2) ) 
IF(Al.LT.O) A I = -A I 
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APPENDIX 4 




DATA 

INPUT FORMAT 

Card 1: 

FORMAT 

(20 A4) 



TITLE 

- Title of the problem 

Card 2: 

FORMAT 

(1015) 


(1) 

INODE 

- No. 

of nodes 

(2) 

NELSM 

- No. 

of elements 

(3) 

NAPC 

- No. 

of applied point 

(4) 

NBC 

- No. 

of free nodes 

Card 3: 

FORMAT 

(1015) 


(1) 

ISW(l) 

= 0, static analysis 


= N, dynamic analysis for N steps 

= -N, static elasto-plastlc analysis for 
N load increments 

(2) ISW(2) = 0, Elastic analysis 


= -1, Viscoelastic analysis 
= 1, Viscoelastoplastic analysis 


(3) 

ISW(3) 

as 

M, Print for each M — time step 

(4) 

ISW(4) 

= 

0 

(5) 

ISW(5) 

= 

No. of rigid elements 

Card 4: 

FORMAT 

(6F13, 

.6) 

(1) 

E 

- 

Modulus of elasticity for rigid element 

(2) 

XNU 

- 

Poisson's ratio for rigid element 

(3) 

DENS 

- 

Soil density 

(4) 

DEPTH 

- 

Maximum soil depth 
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(5) 

ES 

- 

Modulus of elasticity for soil 

(6) 

XNUS 

- 

Poisson's ratio for soil 

Card 5: 

FORMAT 

(6F13 

• 6) 

(1) 

AT 

= 

T( ) x E 0 , where T( r ^ is the soil relaxation 
time in seconds 

(2) 

XK 

- 

Modulus of elasticity for soil 

(3) 

DELT 

- 

Magnitude of time step in seconds 

Ca rd 6 : 

FORMAT 

(6F13, 

.6) 

(1) 

PHI 

- 

Angle of internal friction 

(2) 

VOIDI 

- 

Initial void ratio 

(3) 

CAPA 

- 

Swelling index 

(4) 

RAMD 

- 

Compression index 

Cards 7: 

FORMAT 

(5x , 

2F10.4, 215) 

(1) 

Z(I) 

- 

Z - coordinate value of Node I 

(2) 

R(I) 

- 

R - coordinate value of Node I 

(3) 

IZ 

= 

0 if free to z-direction 



* 

1 if not 

(4) 

IR 

= 

0 if free to R-direction 


$ 1 if not 

Repeat INODE times in the order of node number 
* 

Upward Z is positive 
Cards 8: FORMAT (5x, 515) 

4 corner nodes of an element in counter clockwise. Repeat NELEM 
times in the order of element number. Note that rigid element 
should be numbered first. 

Card(s) 9- FORMAT (I5,2F10.4) 

(1) NODE - Node number with applied load 
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(2) PZ - Z-component 

(3) PR - R-component 

Note that the regular sign convention of theory of elasticity 
is used for input quantify. 


Repeat NAPC times. 
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Initialization 



For incremental elasto-plastic analysis, see App. 3 of Part I 
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APPENDIX 6 

SUBROUTINE ORGANIZATION CHART 




APPENDIX .7 


Subroutine 
Name 

AREAA 

DIAGNL 

DISPL 

F 

FACTOR 

GAUSS 

INIT 

PSTIFF 

SETUP 

SOLTN 

STEP 

STIFF 1 

STIFF 2 

STRAIN 

ZERO 


DESCRIPTIONS OF SUBROUTINES 


Descriptions 


Computes the cross sectional area of an element 

Clears one dimensional array and puts l's on 
diagonal. 

Calls FACTOR and SOLTN, and prints displacement 
vector. 

Functions to be integrated 

Factors (forward substitution) a given simultaneous 
equations 

Integrates by the Gaussian quadrature 

Forms elastic matrix and viscous matrix 

Checks for yielding and forms plastic stiffness 
matrices for yielded elements 

Initializes integration constants 

Backward substitution is performed to give a set of 
solutions to the given simultaneous equations. 

Integrates the equation of motion by step integration 
scheme for dynamic analysis 

Forms elastic stiffness matrix, consistent mass 
matrix, and viscous matrix. Also assembles in 
global form applying the boundary conditions. 

Checks for yielding and forms elasto-plastic stiff- 
ness matrix, and assembles in global form applying 
the boundary conditions. 

Computes strains, stresses, and principal stresses 
(compression is positive). 

Clears any given matrix 


